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ABSTRACT 


A recently proposed taxonomic classification of extant 
ungulates sparked a series of publications that 
criticize the Phylogenetic Species Concept (PSC) 
claiming it to be a particularly poor species concept. 
These opinions reiteratively stated that (1) the two 
fundamental elements of the "PSC", i.e., monophyly 
and diagnosability, do not offer objective criteria as 
to where the line between species should be drawn; 
and (2) that extirpation of populations can lead to 
artificial diagnosability and spurious recognitions of 
species. This sudden eruption of criticism against 
the PSC is misleading. Problems attributed to 
the PSC are common to most approaches and 
concepts that modern systematists employ to establish 
species boundaries. The controversial taxonomic 
propositions that sparked criticism against the PSC are 
indeed highly problematic, not because of the species 
concept upon which they are based, but because no 
evidence (whatsoever) has become public to support 
a substantial portion of the proposed classification. 
We herein discuss these topics using examples 
from mammals. Numerous areas of biological 
research rest upon taxonomic accuracy (including 
conservation biology and biomedical research); hence, 
it is necessary to clarify what are (and what are not) 
the real sources of taxonomic inaccuracy. 


Keywords: Alpha taxonomy; Phylogenetic Species 
Concept; Species concepts; Taxonomic inertia; 
Taxonomic inflation 
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INTRODUCTION 


A recently proposed taxonomic classification for extant 
ungulates (Groves & Grubb, 2011) sparked a series of 
publications criticizing the species concept upon which the 
classification was based, i.e., Phylogenetic Species Concept 
(PSC) (Heller et al., 2013; Zachos et al., 2013; Zachos, 
2013, 2015; Zachos & Lovari, 2013), albeit previous published 
opinions had already presented some of the same arguments 
against the PSC (e.g., Frankham et al., 2012; Garnett & 
Christidis, 2007; Isaac et al., 2004; Tattersall, 2007). Two 
main claims about the PSC have been reiteratively used to 
highlight it as a particularly poor species concept: (1) the 
two fundamental elements of the PSC, i.e., monophyly and 
diagnosability, do not offer objective criteria as to where the 
line between species should be drawn; and (2) the extirpation 
of populations can lead to artificial diagnosability and spurious 
recognitions of species. Moreover, these criticisms portray the 
use of the PSC as detrimental to conservation efforts. We 
argue that the problems attributed to the PSC are common 
to most methodological approaches to species limits and to 
the most commonly used species concepts that have been the 
basis for the taxonomic classifications of mammals currently in 
use. Furthermore, we present evidence that the PSC based 
on diagnosability and monophyly as operational criteria has 
helped to substantially advance mammalian systematics. In 
addition, we show that the recent criticism against Groves & 
Grubb’s (2011) ungulate taxonomy is mistakenly focused on an 
"alleged poverty" of the PSC (see a brief comment relevant to 
this general topic by Tsang et al., 2016, p. 529), whereas the 
real cause of taxonomic inflation in that proposed classification 
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lays on numerous empirical problems. 
PHYLOGENETIC SPECIES CONCEPT (PSC) 


Before we discuss these matters, we must clarify that the 
name “Phylogenetic Species Concept” has been associated 
to various concepts (e.g., McKitrick & Zink, 1988; Nixon 
& Wheeler, 1990), two of which are central in the 
above-mentioned debate. These concepts can be better 
regarded as sets of criteria for species delimitation rather than 
"concepts", as proposed by de Queiroz (2007); however, herein 
we refer to these sets of criteria as "concepts" only to facilitate 
communication by using the same terminology employed by 
authors of previous articles. These concepts are as follows 
(see summaries by Groves et al., 2017; Zachos, 2016a): 


Phylogenetic species concept, diagnosis-based 
version (dPSC): “The smallest diagnosable cluster 
of individual organisms within which there is a 
parental pattern of ancestry and descent" (Cracraft, 
1983, p.170). Subsequently formulated as " 

the smallest aggregation of populations (sexual) 
or lineages (asexual) diagnosable by a unique 
combination of character states in comparable 
individuals (semaphoronts)” (Nixon & Wheeler, 
1990). A more recent version states that species 
are “ the smallest population or aggregation 
of populations which has fixed heritable differences 
from other such populations or aggregations” 
(Groves & Grubb, 2011; see also Groves, 2017). 


Phylogenetic species concept, monophyly-based 
version (mPSC) *... a geographically constrained 
group of individuals with some unique apomorphous 
character, is the unit of evolutionary significance" 
(Rosen, 1978, p. 176). 


A third concept that must be incorporated in the discussion 
is as follows (see Groves et al., 2017): 


Phylogenetic species concept,  diagnosis-and- 
monophyly-based version (dmPSC), defined as 
» the smallest diagnosable cluster of individual 
organisms forming a monophyletic group within which 
there is a parental pattern of ancestry and descent" 
(Mayden, 1997, p. 407; McKitrick & Zink, 1988). 


DO THESE VERSIONS OF THE PSC OFFER OBJECTIVE 
CRITERIA AS TO WHERE THE LINE BETWEEN SPECIES 
SHOULD BE DRAWN? 


The short answer is "no", but "no" would also be the answer 
if the question were asked with regard to any other species 
concept, including the Biological Species Concept (BSC) 
when applied to allopatric populations (see Groves, 2012 and 
references therein). However, the three phylogenetic species 
concepts described above differ importantly regarding the 
degree of objectivity with which they can be applied. 
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With regard to the dPSC (sensu Cracraft, 1983; see also 
Eldredge & Cracraft, 1980; Nixon & Wheeler, 1990; Wheeler 
& Meier, 2000 and references therein), we agree with previous 
criticisms (Heller et al., 2013, 2014; Zachos & Lovari, 2013, 
Zachos, 2016a) in that this concept is prone to promote 
spurious recognition of mere geographic (including subspecies) 
or even individual variants of a single species as if such 
variants were each a valid species. This is due to the high 
degree of subjectivity and arbitrariness implicit in the task 
of judging what characteristics are to be deemed adequate 
to diagnose species and distinguishing such characteristics 
from those that would simply lead to diagnoses of populations, 
or groups thereof, within a single species (but see Wiens 
& Servedio, 2000). Species are not phenotypically and 
genotypically homogeneous across geography, therefore it is 
always the case that various populations within a single species 
can be diagnosed. These diagnoses by themselves must 
not be a justification to regard such populations as different 
species. This limitation of the dPSC is exacerbated when 
sample sizes are small, as is often the case for medium and 
large mammals. In these cases, a researcher may erroneously 
infer the existence of phenotypic discontinuity and the presence 
of characteristics enabling the diagnosis of a sample—the 
latter based on a set of specimens that at the time were 
perceived as worthy of species-level recognition. However, 
as more samples are obtained, individuals with intermediate 
phenotypes with respect to the putative new species and other 
geographic samples may be found. This would render the 
putative new species, which was previously thought to be 
diagnosable, conspecific with an already recognized species 
(e.g., Peres et al., 1996). Examples of these plausible problems 
are abundant in the proposed classification of ungulates by 
Groves & Grubb (2011) (see below), but are by no means 
exclusive to it (e.g., Díaz et al., 1999, 2002; Fonseca & 
Pinto, 2004; Solari, 2004; van Roosmalen et al., 2000, 2007); 
numerous examples exist in early contributions to mammalian 
taxonomy (e.g., Miller, 1912; Pocock, 1941; Robinson & Lyon, 
1901), and even the last decade has seen claims advocating for 
the recognition of a species made on the basis of phenotypic 
diagnoses of as few as one or two specimens—e.g., Meijaard 
et al., 2017 p. 513; see also Mantilla-Meluk (2013) for a 
monkey subspecies named on the basis of morphometric data 
and pelage coloration from only four specimens. Unfortunately, 
in some cases descriptions of species have been carried out 
not only with unacceptably small sample sizes but also merely 
based on images (illustrations, photos, or both) and lacking 
preserved type specimens (see Pine & Gutiérrez, 2018 for a 
review of cases and problems associated to this phenomenon). 
Although no data exist to support the notion that the collection 
of a single individual (for it to properly serve as a preserved 
holotype) significantly increases the probability of an already 
endangered species to become extinct, some researchers may 
prefer not to carry out such collection (e.g., Donegan, 2008; 
but see Dubois & Nemésio, 2007; Dubois, 2009), or it may be 
unfeasible due to impediments in obtaining collection permits. 
In such cases, a wide survey of museum specimens might 


lead to the discovery and subsequent use of specimens in 
taxonomic descriptions. Undertaking comprehensive surveys 
of museum specimens may be disregarded by describers of 
new species, but the possible data yielded, which may be 
coupled with photos of living animals, might ameliorate the 
detrimental effects of extremely small sample sizes and help in 
unveiling geographic and non-geographic (ontogenetic, sexual) 
variation (e.g., Garbino et al., 2016). 

The second concept, the mPSC, under which species must 
be both monophyletic and geographically restricted, seems 
indefensible. Within a species there can be large numbers 
of monophyletic groups that are geographically restricted only 
due to recent changes in their environment. An example 
of this is the populations of brocket deer of the Cordillera 
de Mérida, Venezuela, which for decades were recognized 
as a valid species, Mazama bricenii. This recognition was 
based on no data whatsoever and on the assumption of 
a plausible differentiation due to its supposed geographic 
isolation. However, a recent study that employed ecological 
niche modeling found that, if these populations were truly 
isolated in modern time, such isolation commenced not long 
ago (Gutiérrez et al., 2015). The study also found that while 
the focal population formed a monophyletic haplogroup, it was 
embedded within a larger (yet shallow) clade whose terminal 
branches corresponded to Mazama rufina. Results from that 
study showed that what was once known as Mazama bricenii 
actually corresponds to Mazama rufina (Gutiérrez et al., 2015), 
and illustrate how the application of the mPSC would have led 
to mistakenly recognize M. bricenii as if it were a valid species. 
Such taxonomic recognition would be a mPSC-based artifact 
caused by (1) the fact that sequences obtained from specimens 
from the Cordillera de Mérida (where the populations to which 
the name M. bricenii would apply occur) were recovered in a 
monophyletic haplogroup; and (2) because those populations 
might be geographically isolated in modern time. 

The third concept, dmPSC, for which species must be 
both monophyletic and diagnosable, has been useful for 
improving mammalian taxonomy. As previously noted, many 
monophyletic groups are found within a single species, and 
that monophyly per se does not constitute a criterion to 
determine where the line between species should be drawn. 
However, it is also true that assessing whether a candidate 
species is monophyletic or not provides a fundamental basis 
for its potential recognition as a valid species. Recognizing 
a polyphyletic taxon as if it were a valid species would be 
absurd. On the other hand, in some situations (e.g., species 
originating from peripheral isolation) a candidate species might 
meet the criteria (i.e., monophyly and diagnosability) for validity 
under the dmPSC, but its recognition renders the species in 
which the candidate had thus far been included as paraphyletic. 
No consensus has been reached as to whether taxonomists 
should accept paraphyletic species as valid, or, alternatively, 
if taxonomists should recognize as valid only those that are 
monophyletic (see Carter et al., 2015; Crisp & Chandler, 1996; 
Dias et al., 2005; Ebach et al., 2006; Freudenstein, 1998; 
Funk & Omland, 2003; Hórandl, 2006, 2007; Nelson et al., 


2003; Nixon & Wheeler, 1990; de Queiroz & Donoghue, 1988; 
Zachos, 2014b; Zander, 2007; see also Funk & Omland, 2003). 
Discussing these views requires a much more extensive text 
and would distract from the aim of this perspective piece—i.e., 
clarifying that despite the recent criticisms made against 
PSCs, at least one of these concepts has served to positively 
advance mammalian systematics. By requiring monophyly, 
the application of the dmPSC secures that a phylogenetic 
inference is conducted to describe or revalidate a species, thus 
decreasing the chances that polyphyletic groups of populations 
would be named as a species. These phylogenetic estimates 
also provide frameworks for evaluating alternatives in cases in 
which the description or recognition of a clade as a species 
would render an already recognized species as paraphyletic. In 
such cases, researchers might simply not describe or formally 
recognize that clade at all—which might be acceptable, as 
not every clade in a phylogenetic tree represents a taxon 
worth naming—or describe it at the subspecies level—with 
paraphyly persisting at a lower taxonomic rank—or describe 
it and accept paraphyletic species as valid, in which case the 
researcher could not invoke any species concept that requires 
monophyly (including the dmPSC) upon which to base the 
description. Whichever of these alternatives the researcher 
prefers, and attempts to justify, due to philosophical, pragmatic, 
or both considerations, the fact that the dmPSC requires a 
phylogenetic estimate is an advantage over other concepts that 
do not, including the dPSC and the BSC. 

The requirement that a candidate species must also be 
diagnosable in order to be recognized under the dmPSC is 
indispensable. A species must have a series of genetically fixed 
characteristics that are common to its members and that serve 
to distinguish it from other such species. However, as already 
discussed (see above), diagnosability alone is, in general, an 
inadequate approach to establish species limits. 

Acknowledging the existence of the dmPSC is important 
because it does represent one of the most explicit methods to 
infer species limits—contra authors that ignored the existence 
of this concept in their arguments against or in favor of 
the "PSC" (e.g., Gippoliti et al., 2018; Groves, 2012, 2017; 
Heller et al., 2013; Zachos & Lovari, 2013; Zachos et 
al., 2013; Zachos, 2013, 2014a, 2015, 2016a). Some 
operational steps in delimiting species will always be arbitrary. 
In this sense, the advantage of the dmPSC over other 
concepts is that its operational criteria for recognition of 
species can be objectively tested. In other words, monophyly 
and diagnosability are, in general, more easily testable for 
allopatric populations than reproductive barriers (BSC), and 
more objectively demonstrated than the central criteria upon 
which other concepts define species, such as “ecological 
roles” in the Ecological Species Concept (Van Valen, 1976). 
When applied based on sufficient geographic and taxonomic 
sampling, and, ideally (but not strictly necessary; see 
below), employing phylogenetic inferences using data from 
independent sources (e.g., DNA sequence data obtained from 
independently inherited genes), the dmPSC has improved 
the taxonomic classifications of various groups of mammals, 
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some of which remained problematic for decades. Among 
studies that exemplify how the dmPSC has helped to advance 
mammalian systematics, even if some of them used this 
species concept without explicitly or correctly invoking it, are 
those on didelphid marsupials (e.g., Diaz-Nieto & Voss, 2016; 
Giarla et al., 2010; Gutiérrez et al., 2010; Martinez-Lanfranco 
et al., 2014; Pavan et al., 2017; Voss et al., 2018), rodents (e.g., 
Hawkins et al., 2016; do Prado & Percequillo, 2017; Rogers 
& Gonzalez, 2010; Voss et al., 2013), bats (e.g., Baird et al., 
2008; Molinari et al., 2017; Moras et al., 2016; Velazco et 
al., 2010), and medium and large mammals (e.g., Bornholdt 
et al., 2013; Gutiérrez et al., 2015; Helgen et al., 2009, 2013; 
Janecka et al., 2008; Koepfli et al., 2008; Miranda et al., 
2017; do Nascimento & Feij6, 2017 (and references therein 
for phylogenetic evidence)). These studies have not only 
unraveled the true-species nature of previously unrecognized 
species, but in many cases have shown that taxa considered 
as valid species for decades are not valid species at all. 

The application of any PSC can promote rampant taxonomic 
inflation when applied without sufficient rigor. In our opinion, 
this inflation is caused less by philosophical aspects and 


properties of the PSCs and more by empirical shortcomings. 


In several studies, geographic and individual variation do not 
appear to be satisfactorily addressed, and names are applied to 
what could be intraspecific variants (see examples in Tattersall, 
2007). On other occasions, monophyletic groups recovered from 
molecular phylogenies based on sequence data from a single 
locus promptly receive new or revalidated names (e.g., Boubli 
et al., 2012; Thinh et al., 2010). Nevertheless, it is important to 
note that when the established, traditional taxonomic classification 
of a focal group is the result of dogmatic acceptance of expert 
opinions (often past-century authorities), without support from 
data (see Gutiérrez & Helgen, 2013), then even the use of 
limited evidence—e.g., analyses of sequence data from a single 
locus (despite the well-known shortcomings of this approach; see 
Dávalos & Russell, 2014; Knowles & Carstens, 2007; Maddison, 
1997), ideally coupled with qualitative and/or quantitative analyses 
of morphological data—can well justify taxonomic changes if 
based on adequate sampling (e.g., Gutiérrez et al., 2010, 2015, 
2017; Voss et al., 2013; contra Zachos, 2009, 2016b). 

The proposed ungulate taxonomic classification that 
sparked the recent series of criticism against the PSC 
is particularly problematic, but not because it was based 
on the dPSC. Most controversial aspects of this proposed 
classification are not at all associated to any species concept, 
as it might seem if one reads the recent debate between Frank 
Zachos and Colin Groves and their co-authors with regard 
to the "PSC", but rather to more practical aspects of such a 
monograph, to name just a few (see also Heller et al., 2013; 
Holbrook, 2013; Zachos, 2014a, p. 1): (1) Groves & Grubb 
(2011) did not assess geographic variation at all for most of the 
species they recognized; (2) unfortunately, for some species 
recognized by these authors, the sample size employed was 
not indicated nor any published study cited to support the 
taxonomic proposals, whereas for many other alleged species 
the sample sizes were extremely low—e.g., Alcelaphus tora, 
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Dorcatragus megalotis, Eudorcas nasalis, Eudorcas tilonura, 
Gazella acaciae, Gazella karamii, Gazella shikarii, Lama 
mensalis, Madoqua hararensis, Mazama fuscata, Mazama 
jucunda, Mazama trinitatis, and Redunca cottoni; (3) no 
published phylogenetic information seems to be the basis of 
most of their taxonomic propositions; (4) in general, no detailed 
discussions were presented on whether recognizing a taxon 
as a valid species was more appropriate and justifiable than 
regarding it as a subspecies, and it seems that the objective 
of the authors was to merely recognize as valid species as 
many taxa as possible, without critical evaluation of alternatives 
(e.g., recognizing subspecies when appropriate); (5) a list 
of the specimens examined was not provided, and hence it 
is difficult for the scientific community to evaluate the authors' 
assertions on specimen morphologies based on the same material 
with certainty; (6) no data were made available that would enable 
reproduction and testability of the analyses that were the basis of 
the taxonomic propositions; (7) although the authors cite published 
studies for some of the taxonomic changes they proposed, for 
others they did not and nowhere in their monograph can be found 
results from any quantitative analyses. We cannot understand 
why Groves & Grubb (2011) failed to publish the results of their 
quantitative analyses given current possibilities to do so (see below). 
Unfortunately, this problem is not unique to the proposed ungulate 
taxonomic classification of Groves & Grubb. An important volume 
on mammals of South America (Patton et al., 2015) contained 
the first modern taxonomic treatments of various rodent groups 
(e.g., family Sciuridae, genera Aepeomys, Oecomys, Rhipidomys, 
Thomasomys), but results from analytical procedures assessing 
geographic and non-geographic variation of those groups have not 
been published (in the book or elsewhere), and in some cases 
it seems unlikely they will ever be published. Luckily, several 
unpublished Ph.D. dissertations that served as the basis for the 
book sections treating those taxa have been privately shared 
among colleagues. Clearly, making these Ph.D. dissertations 
(and other unpublished material) digitally available to the scientific 
community free of charge from a repository on the Internet (e.g., 
Dryad, Figshare, Internet Archive, ResearchGate, Zenodo), if their 
authors grant authorization, should be considered by the editors of 
this book, and similar actions should be considered by authors and 
editors of future monographs introducing taxonomic classifications. 


CAN EXTIRPATION OF POPULATIONS LEAD TO ARTIFICIAL 
DIAGNOSABILITY AND SPURIOUS RECOGNITIONS OF 
SPECIES UNDER THE dmPSC? 


The short answer is "it can", but again the same answer 
would apply if the question were asked about other species 
concepts, including the BSC. In their criticism of the "PSC", 
Zachos & Lovari (2013) claimed that the fact that extirpation of 
populations can lead to artifactual diagnosability and monophyly 
is one of the weaknesses of the PSC that makes this concept 
a particularly poor one. They stated that “There is yet another 
line of argumentation that clearly shows the shortcomings of 
both diagnosability and monophyly as yardsticks for species 
delimitation and that we believe is another coup de grace for the 
PSC. Diagnosability (just like reciprocal monophyly) can and 


often does occur as a consequence of extinction of intermediate 
forms...”. Unfortunately, Zachos & Lovari (2013) did not realize 
that extirpation of intermediate populations is one of the natural 
causes of speciation. These events lead to speciation, affecting 
gene flow and ecological adaptation of extreme phenotypes in 
populations that previously were genetically connected by the 
existence of intermediate populations. In fact, it could be said 
that living taxa that can be validly recognized on the planet exist 
as separate biological entities and as taxonomically diagnosable 
units only because of the extinction of intermediate forms. If 
all intermediate forms that have lived since the beginning of 
life on Earth were still with us, then all living organisms, from 
prokaryotes to eukaryotes and from plants to animals, would 
exist as a single morphological and reproductive continuum: 
no distinguishable taxa would exist! Logically, in cases in 
which extirpations have taken place fairly recently (e.g., due to 
human-related causes) no speciation may have yet occurred. 
In such cases, the extirpation of populations can potentially 
lead to artifactual recognitions of species under the PSCs, but 
this issue is far from being associated only to the PSCs; rather, 
it is a problem that can potentially affect most, if not all, 
species concepts currently in use. For example, this issue 
can lead to artifactual recognition of species under the BSC. 


Let us imagine a species, species A, with wide distribution 
and showing geographic variation by way of a cline in several 
qualitative cranial traits believed to be taxonomically important, 
i.e., used by most authors to distinguish species within the 
corresponding genus. We will illustrate these traits as color 
and shape in polygons that represent populations of species 
A (Figure 1, panel 1). If recent extirpations of intermediate 
populations take place and only those populations occurring 
at the opposite extremes of species A’s range remain as extant, 
then these populations will become allopatric and it would be 
highly likely that they would be considered as members of 
different species under the BSC (Figure 1, panel 2). The BSC 
would fail to regard these populations as conspecific, even 
employing the approach presented by Tobias et al. (2010; 
see also Brooks & Helgen, 2010), which uses the degree of 
differentiation known to exist between different but sympatric 
species (i.e., species A and B in Figure 1) as a standard 
to assess the taxonomic relevance of differentiation between 
allopatric populations (i.e., populations of species A in North 
and South America in Figure 1)—in order words, this approach 
uses the former degree of differentiation as a threshold at which 
(or above) allopatric populations could be treated as different 
species under the BSC. 





Figure 1 Illustration of how population extirpations can promote artifactual recognition of populations of a single species as if they 
were different species under the Biological Species Concept (and many other concepts) 


Species A, which is represented by quadrilaterals (rectangle and rectangle-like polygons), possesses a wide distribution and geographic variation by way of a cline in 


several qualitative cranial traits considered taxonomically important; herein these traits are illustrated as color and shape of the polygons (panel 1). Species B, which 


is represented by circles, is restricted to North America, where it occurs in sympatry with some North American populations of species A (panel 1). If extirpation of 


populations of species A takes place and only those populations occurring at the opposite extremes of species A's range remain as extant, then these populations will 


become allopatric (panel 2). In that scenario, it would be highly likely that the remaining extant populations would be considered as members of a different species 


under the BSC (panel 2). The BSC would fail to regard these populations as conspecific, even employing the approach presented by Tobias et al. (2010; see also 


Brooks & Helgen, 2010), which uses the degree of differentiation known to exist between different but sympatric species (i.e., species A and B) as a standard to 


assess the taxonomic relevance of differentiation between allopatric populations (i.e., populations of species A in North and South America). 


CONCLUSIONS 


Major differences exist among the different concepts labeled 
as the "Phylogenetic Species Concept", and the one that 


uses both diagnosis and monophyly (dmPSC) to delimit 
species has been, and will continue to be, important for 
positively advancing mammalian taxonomy. Our preceding 
discussion should rectify misunderstandings that could arise 


Zoological Research 39(5): 301—308, 2018 305 


from claims made in recently published opinions debating 
alleged pros and cons of the "PSC". Although we partially 
agree with some of the arguments presented by participants 
of that debate, the proposed taxonomic classification of 
ungulates (Groves & Grubb, 2011) that motivated this debate 
is highly deficient, in our view, not so much because of the 
species concept it employed (i.e., dPSC), but rather due to 
serious empirical problems. Among them is the absence 
of statistical assessments of geographic and non-geographic 
variation in diagnostic traits. Although in many instances 
the number of specimens available in museums should have 
permitted statistically satisfactory assessments of geographic 
and non-geographic variation, results from those analyses 
were not presented by Groves & Grubb (2011) in their 
monograph. In other instances, extremely low sample sizes 
precluded proper statistical analyses. 

To refrain from producing taxonomic hypotheses because 
of limited material (e.g., few available museum specimens) 
would hamper progress in medium and large mammal 
taxonomy. As previously mentioned, collecting new samples 
of such mammals can be logistically impracticable, and some 
researchers may simply prefer not to collect them due to 
conservation concerns (but see clarification above). Thus, 
even when the available material consists of only a few 
specimens, taxonomic studies should still be carried out, but 
the taxonomist should bear in mind the statistical limitations of 
a small sample, such as inadequate estimations of population 
ranges and lower confidence levels. Collating information 
from multiple data sources, such as nucleotide sequences, 
discrete and continuous morphological data, and behavior—an 
approach nowadays called “integrative taxonomy” (e.g., Dayrat, 
2005)—has long been considered useful as it theoretically 
increases the probability of correctly identifying and delimiting 
taxonomic entities (Simpson, 1961; Tinbergen, 1959). 

Numerous areas of biological research rest upon taxonomic 
accuracy (including conservation biology and biomedical 
research); hence, it is necessary to clarify what are (and what 
are not) the real sources of taxonomic inaccuracy. 
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ABSTRACT 


Apodemus (mice) and Rattus (rats) are the top rodent 
reservoirs for zoonoses in China, yet little is known 
about their diversity. We reexamined the alpha 
diversity of these two genera based on a new collection 
of specimens from China and their cyt b sequences 
in GenBank. We also tested whether species 
could be identified using external and craniodental 
measurements exclusively. Measurements from 147 
specimens of Apodemus and 233 specimens of 
Rattus were used for morphological comparisons. We 
analysed 74 cyt b sequences of Apodemus and 100 
cyt b sequences of Rattus to facilitate phylogenetic 
estimations. Results demonstrated that nine species 
of Apodemus and seven species of Hattus, plus a 
new subspecies of Hattus nitidus, are distributed 
in China. Principal component analysis using 
external and craniodental measurements revealed 
that measurements alone could not separate the 
recognized species. The occurrence of Rattus pyctoris 
in China remains uncertain. 


Keywords: Alpha diversity; Apodemus; DNA-barcoding; 
Rattus; Taxonomy; Phylogenies; New subspecies 


INTRODUCTION 


Small volant and nonvolant mammals are important 
components of ecological communities and play vital roles in 
ecological systems. They are among the most common agents 


for infections and, thus, have strongly affected human history. 


For example, black rats (Rattus rattus) are considered likely 
agents for the spread of Oriental rat fleas, which drove the 
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Black Death plague throughout Europe and the Mediterranean 
during the 14th century and killed 30%—60% of the European 
population (Barnett, 2001; Duplantier et al., 2003). More 
recent examples of small mammal zoonoses include severe 
acute respiratory syndrome (SARS) caused by a coronavirus 
and Ebola hemorrhagic fever caused by Ebolavirus, with hosts 
including, but not limited to, bats and civets (Klein & Calisher, 
2007; Menachery et al., 2015). Rodent-borne diseases 
such as plague and hantavirus have made considerable 
contributions to human illnesses and are responsible for more 
deaths than all wars combined (Klein & Calisher, 2007). New 
pathogens, especially hantaviruses, have been isolated from 
rodents in China and adjacent countries annually (Huang et 
al., 2017). Because different species have specific immune 
systems and different levels of tolerance to zoonotic infections, 
identification of rodent reservoirs of zoonotic pathogens is a 
high priority (Meerburg et al., 2009). 

Rats and mice often top the zoonoses reservoir list of the 
Chinese Center for Disease Control and Prevention (China 
CDO) because of the large number of species, substantial 
population sizes, and high potential for carrying zoonotic 
pathogens (Wu et al, 2017). Unfortunately, we still do 
not know how many species of rats and mice occur in 
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China, or which species carry what pathogens, even for the 
most common genera such as Apodemus and Rattus. The 
reasons for this are complicated. Both Apodemus and Rattus 
have complex evolutionary and taxonomic histories, with 
classifications continuously being updated. Switching between 
valid species and synonyms causes considerable confusion, 
especially for non-specialist researchers. Furthermore, many 
species occur only in remote mountains or near national 
borders with high species diversity, such as Yunnan, Xizang 
(Tibet), and Xinjiang. Indeed, the rats and mice of southern 


Xizang and western Xinjiang remain to be studied carefully. 


Finally, many rodents are difficult to identify to species level 
due to the number of morphologically similar species (Galan et 
al., 2012). 

The latest version of Mammal Species of the World (Musser 
& Carleton, 2005) recognized 20 species of Apodemus and 
162 synonyms. Several scenarios for the classification of 
Apodemus have been proposed (Filippucci, 1992; Martin et al., 
2000; Musser et al., 1995; Serizawa et al., 2000; Zimmermann, 
1962), but none are strongly supported, and phylogenies 
remain poorly resolved despite molecular efforts (Liu et al., 
2004; Serizawa et al., 2000; Suzuki et al., 2003). Furthermore, 
the number of species in China remains unknown, with 
previous estimations varying from six (Corbet, 1978; Xia, 1984), 
seven (Liu et al., 2002; Liu et al., 2004), eight (Smith et al., 
2008), and nine species (Nowak, 1999; Wang, 2003). Many 
authors have suggested that A. sylvaticus occurs in Xinjiang, 
China (Corbet, 1978; Xia, 1984; Wang, 2003), whereas others 


have argued that the species is A. uralensis (Smith et al., 2008). 


The former species occurs in Western Europe (Bousbouras, 
1999; Macholán et al., 2001; Mezhzherin & Zykov, 1991; 
Michaux et al., 1996), and its incorrect identification in China 
likely relates to outdated taxonomy. 

Rattus, another problematic genus, has had 25 subgenera and 


more than 550 species and subspecies named (Simpson, 1945). 
Currently, 66 species are recognized but uncertainty persists. 


Previous supermatrix analysis did not obtain a monophyletic 
Rattus, indicating that systematics is far from resolved (Steppan 
& Schenk, 2017). Arguments also persist for the most common 
species, including black rats whose species boundary remains 
unfixed (Aplin et al., 2011). The number of species of Rattus 
in China is also uncertain and varies from four (Corbet, 1978), 
seven (Smith et al., 2008), and nine (Wang, 2003). 

Similar to other rodents, species in these two genera 
are difficult to identify or distinguish morphologically due to 
their similar appearance, overlapping measurements, and key 
factors involving the single cusp on their teeth. Diagnosis 
often requires clean skulls, which are not always available or 
correctly prepared. DNA barcoding is a promising approach 
but requires a solid reference database (Moritz & Cicero, 
2004). Unfortunately, GenBank data are problematic because 
many rodent sequences are uploaded by non-specialists 
such as epidemiological researchers. This reduces the 
reliability of environmental assessment reports and hampers 
our understanding of host and disease associations. 

Herein, we revisited the alpha diversity of Apodemus and 
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Rattus in China based on a collection of more than 400 
specimens and the integration of cyt b sequences. We 
evaluated the species of both genera in China and assessed 
if they could be identified easily using traditional morphometric 
approaches. 


MATERIALS AND METHODS 


Morphological diagnoses and analyses 

We examined 147 specimens of Apodemus and 233 
specimens of Rattus collected from multiple localities across 
China. External and skull measurements followed Liu et 
al. (2012). External measurements of fresh specimens in 
the field were taken to the nearest 0.5 mm using a steel 
tape. These included head-body length (HBL), hind-foot length 
(HFL), ear length (EL) and tail length (TL) (museum specimens 
from original records). We measured eight skull variables 
using a digital caliper graduated to the nearest 0.01 mm 
from 147 intact skulls of Apodemus and 233 intact skulls 
of Rattus, including greatest length of skull (SGL), nasal 
bone length (NBL), zygomatic breadth (ZB), skull basal length 
(SBL), upper toothrow length (UTRL), upper molar row length 
(UMRL), auditory bulla length (ABL), and mandible length 
(ML). Examined specimens (Supplementary Table S1) were 
deposited in the Kunming Institute of Zoology (KIZ), Sichuan 
Academy of Forestry (SAF), Beijing Institute of Zoology (BIZ), 
Guangdong Key Laboratory of Animal Conservation and 
Resource Utilization, and Fujian Center for Disease Control 
Prevention. 

Specimens were roughly identified based on external and 
craniodental morphology, following Kaneko (2010) and Smith 
et al. (2008). External and craniodental measurements largely 
overlapped between species (see Results) and were inadequate 
for identification. However, several diagnostic characters on 
the upper molars were constructive in classification, including 
the number of lingual angles of the first and second upper 
molar, presence/absence of cusp t3 on the first upper molar, 
and numbers of internal lobes on the third upper molars. We 
also cross-checked results based on morphological diagnoses 
with molecular sequences (when available) to refine identification. 
All specimens were identified by the same researcher (SYL) for 
consistency. We finally assigned our specimens to nine species 
of Apodemus, seven species of Rattus, and a new subspecies of 
Rattus nitidus, respectively. 

We analyzed morphometric variation using principal 
component analyses (PCAs) on log10-transformed variables 
using two datasets for each species. The first dataset included 
both external and craniomandibular variables, whereas 
the second dataset included craniomandibular variables 
only. Inclusion of the external data tested whether these 
measurements could increase the accuracy of identification. 
Statistical analysis was performed using SPSS v16.0 (SPSS 
Inc., USA). When two or more recognized species were not 
well separated in the principal component (PC) plots, analysis 
of variance (ANOVA) was applied to analyze among group 
differences. 


Molecular analysis 

We sequenced mitochondrial cyt b for 74 and 100 
specimens of Apodemus and Hattus, respectively. Localities 
of molecular samples used from China are mapped in 
Figure 1. All sequenced specimens were deposited in the 
SAF. Total genomic DNA was extracted using the standard 
phenol-chloroform method (Sambrook & Russell, 2001). We 
used the universal primers of mammalian cyt b L14724/H15915 
for amplification (Irwin et al., 1991). Polymerase chain reaction 
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(PCR) was conducted in a 25-uL reaction volume, including 2.5 
uL of 10xEX Taq buffer (Mg?* Free), 2 uL of 2.5 mmol/L dNTP, 
1.5 uL of 25 mmol/L MgCls, 1 uL of 10 umol/L primers, and 
1 unit of EX Taq polymerase (TaKaRa Biotech, Dalian, China). 
The product was purified using an EZNA™ Gel Extraction Kit 
(Omega, USA), and was sequenced using the same primers for 
amplification on an ABI 3730XL sequencer. Sequences were 
assembled and edited using SeqMan and EditSeq (DNASTAR, 
Lasergene v7.1) before subsequent analyses. 
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Figure 1 Localities of molecular samples from China in this study 


To avoid misidentification, we first conducted a "naive 
identification" for the obtained sequences using the "identify 
organism" workflow in Geneious v11 (Biomatters, New 
Zealand). The software blasted each sequence against the 
GenBank nucleotide collection (nr/nt) database. When pairwise 
identity between the query (our sequence) and subject (in 
GenBank) sequences was higher 9895, Geneious considered 
them as the same species. We cross-checked the results of 
both morphological and molecular identifications, and when the 
identification was inconsistent, we revisited the skin and skull 
specimens before applying an identity. 

To provide a better picture of species diversity in China, we 
downloaded cyt b sequences of Apodemus (n=477) and Rattus 
(n2273) in China from GenBank, discarding sequences shorter 
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than 800 bp. We also included cyt b data representing another 
12 species of Apodemus and 14 species of Rattus from outside 
of China. An additional five sequences of H. pyctoris from 
Nepal were included. For better estimation of phylogenetic 
relationships, we downloaded the mitochondrial genomes 
(mitogenomes) of seven species of Apodemus and 14 species 
of Rattus (Supplementary Table S2). One mitogenome under 
the name of "Apodemus chejuensis" may not have been a valid 
species. Cyt b of Tokudaia spp. (n=3) and a mitogenome of 
Bandicota indica were selected as outgroup representatives 
for Apodemus and Rattus, respectively, following Steppan & 
Schenk (2017). In total, the datasets for Apodemus and Rattus 
included 572 (with seven mitogenomes) and 397 sequences 
(15 mitogenomes), respectively. We aligned the sequences for 
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each genus using MAFFT v7.3 implemented in Geneious v11. 
We removed all tRNAs, D-loop, and ND6 sequences from the 
alignments, and only used rRNAs and 13 protein-coding genes 
for phylogenetic analyses. Sequence genetic distances were 
calculated for cyt b using MEGA v.5 (Tamura et al., 2011) under 
the Kimura 2-parameter model (Kimura, 1980). 


Phylogenetic analyses 

We employed RAxML v8.2.10, a maximum likelihood-based 
approach, for phylogenetic analyses. We partitioned the 
alignments by genes, except for cyt b, which we partitioned into 
the 15t+21d ang 3 codon positions. Analyses were performed 
on the CIPRES Science Gateway. We used GTR+G as the 
evolutionary model for each partition because RAxML does 
not accept models other than GTR or GTR+G. We ran each 
analysis using the rapid bootstrapping algorithm and let RAxML 
halt bootstrapping automatically. We also repeated analyses 
using alternative strategies, such as different partitioning 
schemes (e.g., partitioned by gene and codon positions 
for all coding genes) and evolutionary models (e.g., using 
GTR model instead of GTR+G), none of which strongly 
altered phylogenetic relationships (i.e., different relationships 
supported by bootstrap values (BS) >75)). 


RESULTS 


Morphological analysis 

Morphological analysis of Apodemus 

Morphological measurement statistics of the eight Apodemus 
species, excluding A. semotus, are given in Table 1. In 
the first PCA, using all 12 measurements (n=139), the 
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first and second principal components accounted for 57.6% 
(eigenvalue=6.9; Table 2, a) and 11.7% (eigenvalue=1.4) of 
total variation, respectively, with all other principal components 
having eigenvalues smaller than 1. PC1 was positively 
correlated with all craniodental variables (loadings>0.63), and 
PC2 was positively correlated with external measurements 
(loadings>0.55). The PC1 and PC2 plot (Figure 2A) did not 
clearly separate the species. Apodemus latronum plotted on 
the positive regions of PC1 and PC2, indicating a large body, 
long tail, long hindfeet, and long ears. In accordance with 
its small skull and small external measurements, A. uralensis 
occurred along the negative regions of PC1 and PC2. The 
sister- or closely related species A. agrarius and A. chevrieri 
as well as A. pallipes and A. uralensis were well separated, but 
both pairs overlapped with A. peninsulae, A. draco, and A. ilex, 
which, in turn, largely overlapped. For the second PCA, using 
eight craniodental measurements (n=141), the first principal 
component accounted for 69.5% of variation (eigenvalue=5.6; 
Table 2, b). The other principal components accounted for 
less than 9.4% (eigenvalue« 0.75) of total variation, indicating 
they were not stable (Shankardass, 2000). Seven variables 
were positively correlated with PC1 (loading>0.56), except for 
UMRL (loading=0.076), which was positively correlated with 
PC2 (loading=0.93). The PC1 and PC2 plot (Figure 2B) was 
similar to the previous plot. None of these species were clearly 
separated from all others. Apodemus chevrieri and A. latronum 
plotted on the positive regions of PC1, indicating a relatively 
large skull, and A. uralensis occurred along the negative region 
of PC1 in accordance with its small skull. 





Wm o e 
0 =. " 4 e 7. 6 
vc e 
O e * "wo. ". 
à e e e 
e e o 
e e 
* — * e 
* 
* * 
* *+ A 
-2.0 F 
e 
+ 
@ A. agrarius 
* 9 A. chevrieri 
W A. draco 
E A. ilex 
B A. latronum 
" * A. pallipes 

9 A. peninsulae 

-4.0 H ** A. uralensis 

















-3.0 -2.0 -1.0 0 1.0 20 3.0 


Figure 2 Principal component analysis of the first three principal components among the eight species of Apodemus based on 


both external and craniomandibular variables (A) and craniomandibular variables only (B) 
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Table 1 Measurement statistics of Apodemus 
SGL NBL ZB SBL 





n Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum 























A. latronum 15 28.03 2.54 20.23 30.86 11.12 0.60 9.70 12.17 13.41 0.67 12.38 14.96 25.21 1.58 22.37 27.52 
A. peninsulae 38 27.76 1.60 24.31 31.75 10.64 1.02 8.09 12.42 13.49 0.96 12.01 15.35 25.58 1.87 22.50 29.81 
A. chevrieri 8 29.02 0.83 27.83 29.99 10.65 0.74 9.16 11.63 13.63 0.34 13.06 14.02 26.29 1.24 24.47 27.57 
A. agrarius 41 26.57 0.97 24.55 28.51 9.55 0.53 8.50 10.67 12.56 0.52 11.23 13.34 24.65 1.20 20.82 26.73 
A. ilex 10 26.63 1.25 25.31 28.88 10.10 0.79 9.13 11.44 12.61 0.37 12.09 13.31 23.42 1.19 22.17 25.45 
A. draco 9 26.80 1.63 23.83 29.19 10.56 1.06 8.65 11.80 12.59 0.52 11.83 13.27 23.64 1.59 20.97 25.41 
A. pallipes 12 26.84 0.58 25.76 27.98 10.37 0.37 9.45 10.95 12.73 0.26 12.26 13.13 23.63 0.56 22.37 24.30 
A. uralensis 14 2437 0.89 22.49 25.45 8.72 0.48 7.96 9.74 12.38 0.52 11.68 13.20 22.66 0.96 20.49 23.82 
UTRL UOSL TBL ML 
n Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum 
A. latronum 15 14.41 0.76 12.94 15.60 5.00 0.23 4.47 5.39 5.52 0.35 4.90 6.38 10.48 3.58 6.16 14.41 
A. peninsulae 38 13.71 0.94 11.52 15.82 4.27 0.29 3.41 5.10 5.34 0.40 4.53 6.17 13.25 1.11 11.02 15.79 
A. chevrieri 8 14.55 0.64 13.57 15.52 4.59 0.23 4.24 4.91 5.90 0.25 5.52 6.19 14.61 0.67 13.70 15.74 
A. agrarius 41 12.97 0.73 10.99 14.64 4.18 0.27 3.45 4.84 5.22 0.31 4.55 5.69 12.61 0.70 10.28 14.06 
A. ilex 10 12.74 0.66 11.98 13.95 4.42 0.37 3.97 4.97 5.33 0.25 4.90 5.73 13.21 0.83 11.97 14.97 
A. draco 9 13.10 0.72 11.98 14.03 4.24 0.06 4.15 4.35 5.23 0.35 4.50 5.62 13.03 1.09 11.36 14.32 
A. pallipes 12 13.25 0.32 12.73 13.83 4.34 0.16 4.14 4.67 4.83 0.13 4.62 5.12 13.51 0.40 12.98 14.10 
A. uralensis 14 11.23 0.47 10.51 12.00 3.60 0.25 2.89 3.82 4.68 0.29 4.30 5.23 11.57 0.50 10.50 12.28 
HBL TL HFL EL 
n Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum 
A. latronum 15 100.90 9.84 88.00 115.00 104.27 8.41 90.00 122.00 24.17 0.98 22.50 26.00 20.20 0.80 19.00 22.00 
A. peninsulae 38 100.32 11.96 73.00 137.00 94.26 10.96 72.00 121.00 22.95 1.59 18.00 25.00 16.08 1.58 14.00 21.00 
A. chevrieri 8 102.50 8.78 90.00 113.00 94.12 9.40 82.00 105.00 23.63 1.69 21.00 26.00 16.12 1.89 13.00 18.00 
A. agrarius 41 100.17 8.28 85.00 120.00 82.54 11.55 50.00 101.00 19.55 1.86 14.00 24.00 13.39 1.05 12.00 16.00 
A. ilex 10 86.30 5.83 80.00 98.00 92.80 11.24 72.00 107.00 21.90 1.10 20.00 24.00 17.50 1.43 15.00 19.00 
A. draco 9 92.11 9.85 75.00 103.00 105.11 12.08 84.00 119.00 48.89 77.30 21.00 255.00 17.33 0.71 16.00 18.00 
A. pallipes 12 93.75 3.93 87.00 99.00 95.67 4.14 90.00 104.00 21.58 0.67 21.00 23.00 16.58 0.52 16.00 17.00 
A. uralensis 14 83.86 7.31 67.00 94.00 74.36 4.99 65.00 83.00 19.82 1.14 18.00 22.00 12.11 2.65 9.00 17.00 


Abbreviations are explained in the Materials and Methods section. All measurements are in mm. 
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One-way analysis of variance (ANOVA) revealed that the seven 
species differed significantly (P<0.05) in all external and cranial 
characters tested, except for NBL (P=0.497), TL (P=0.064), and 
HFL (P=0.094). Results showed significant differences as follows: 
UTRL, MRL, ABL, and ML between A. peninsulae and A. chevrieri; 
ZB, SBL, UTRL, HBL, and EL between A. peninsulae and A. ilex; 
SGL, ZB, SBL, HBL, TL, HFL, and EL between A. peninsulae 
and A. draco; ZB, SBL, and ABL between A. peninsulae and A. 
pallipes; SGL, ZB, SBL, UTRL, ABL, ML, HBL, and EL between A. 
chevrieri and A. ilex; SGL, ZB, SBL, MRL, UTRL, ABL, ML, HBL, 
TL, and HFL between A. chevrieri and A. draco; SGL, ZB, SBL, 
UTRL, MRL, ABL, and ML between A. chevrieri and A. pallipes; 
TL and HFL between A. ilex and A. draco; and ABL, TL, and HFL 
between A. draco and A. pallipes. Thus, morphological analysis 
indicated that the eight species of Apodemus could be separated 
by the 12 morphological characters, validating the taxonomic 
status of these species in China. 


Table 2 Factor loadings and percentage of variance explained 
for principal component analysis 














Apodemus Rattus 
a b c d 
Variables PCI PC2 PC1 PC2 PCi PC2 PCi  PC2 
SGL 0.82 0.07 0.83 0.08 0.96 -0.05 0.97 -0.03 
NBL 0.64 0.39 0.84 0.08 089 -0.06 0.90 -0.03 
ZB 0.83 0.00 0.99 -0.24 0.92 -0.08 0.92 -0.08 
SBL 1.01 -0.17 0.95 -0.09 0.97 -0.08 0.97 0.01 
UTRL 0.80 0.23 0.72 037 0.95 -0.16 0.97 -0.19 
UMRL 029 055 0.08 093 O71 0.15 0.71 -0.35 
ABL 0.77 0.00 0.56 0.34 068 0.14 0.70 0.65 
ML 0.64 024 0.71 017 092 -0.14 0.93 0.02 
HBL 0.09 -0.09 NA NA 0.79 -0.25 N/A N/A 
TL 013 074 NA NA 067 052 NA NA 
HZL -0.04 0.82 N/A N/A 0.74 -0.24 N/A N/A 
EL -0.09 0.94 N/A N/A 0.40 084 NA N/A 
Eigenvalues 6.9 1.4 5.60 0.75 799 1.21 632 057 
Total variance ES 414.7 69.50 9.40 66.54 10.08 7905 7.08 


explained (%) 





For abbreviations see Materials and Methods. N/A: Not available. 


Morphological analysis of Rattus 

Morphological measurement statistics of the seven species of 
Rattus and a putatively new subspecies of R. nitidus (from 
southern Xizang) are given in Table 3. In the first PCA, which 
used all 12 measurements (n=233), the first and second principal 
components accounted for 66.54% (eigenvalue=7.99) and 10.08% 
(eigenvalue=1.21) of total variation, respectively (Table 2, c), with 
all other principal components having eigenvalues smaller than 
1. Most species largely overlapped (Figure 3A). In the second 
PCA, which used eight craniodental measurements (n-233), the 
first and second principal components accounted for 79.05% 
(eigenvalue=6.32) and 7.08% (eigenvalue=0.57) of total variation, 
respectively (Table 2, d). The PC1 and PC2 plot (Figure 3B) 
revealed largely overlapping species. 

One-way ANOVA demonstrated significant differences 
(P<0.05) between the seven species in all external and cranial 
characters tested. Results showed significant differences as 
follows: SGL, NBL, ZB, SBL, UTRL, UMRL, ML, BL, HBL, HFL, 
and EL between H. nitidus and R. losea; NBL, UTRL, and HFL 
between R. nitidus and R. tanezumi, HBL and HFL between 
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H. nitidus and R. andamanensis; BL, HBL, and EL between H. 
nitidus and R. norvegicus; SGL, NBL, ZB, SBL, UTRL, UMRL, 
ML, BL, HFL, and EL between R. nitidus and R. exulans; SGL, 
ZB, SBL, UMRL, HBL, HFL, and EL between R. /osea and R. 
tanezumi; SGL, NBL, ZB, SBL, UTRL, UMRL, ML, BL, HBL, 
HFL, and EL between H. /osea and R. andamanensis; SGL, 
NBL, ZB, SBL, UTRL, UMRL, ML, BL, and HFL between H. 
losea and R. norvegicus; ML, BL, and HFL between H. losea 
and R. exulans; in HBL, HFL, and EL between R. /osea and 
H. rattus; UTRL between H. tanezumi and R. andamanensis; 
UTRL, ML, BL, HBL, HFL, and EL between R. tanezumi and 
R. norvegicus; SGL, NBL, ZB, SBL, UTRL, UMRL, ML, BL, 
HBL, HFL, and EL between H. tanezumi and R. exulans; HBL, 
HFL, and EL between H. andamanensis and H. norvegicus; 
SGL, NBL, ZB, SBL, UTRL, UMRL, ML, BL, HBL, HFL, and EL 
between R. andamanensis and R. exulans; UMRL between R. 
andamanensis and R. rattus; SGL, NBL, ZB, SBL, UTRL, UMRL, 
ABL, ML, and HFL between R. norvegicus and R. exulans; HBL 
and EL between R. norvegicus and R. rattus; and SGL, SBL, 
BL, HBL, HFL, and EL between R. exulans and R. rattus. Thus, 
the 12 morphological characters separated the seven species of 
Rattus and validated their occurrence in China. 

When all individuals of the two subspecies of R. nitidus were 
subjected to an independent sample t-test for each variable, 
significant differences appeared in UTRL, UMRL, ML, and TL 
between R. nitidus nitidus and R. nitidus from Xizang. 


Molecular analysis 

We obtained cyt b sequences for 78 specimens of 
Apodemus and 106 specimens of Rattus. De novo 
sequences were deposited in GenBank under accession Nos. 
MG748165-MG748348 (Supplementary Table S3). 

Cyt b K2P interspecies distances for Apodemus ranged 
from 5.4% to 20.7% (Supplementary Table S4). The smallest 
distance occurred between A. uralensis and A. pallipes, and 
largest between A. sylvaticus and A. latronum. The distances 
for Rattus ranged from 2.196 to 16.5% (Supplementary Table 
S5). The smallest distance occurred between H. baluensis 
and R. tiomanicus, and the largest between R. leucopus and 
R. argentiventer. The K2P distance of R. nitidus from Xizang 
and H. nitidus nitidus was 0.019. 


Matrilineal genealogy (haplotype phylogeny) of Apodemus 
Matrilineal genealogy using the mitogenome and cyt b 
data for Apodemus (n-569) did not fully resolve the 
higher relationships (Figure 4), as in previous studies (see 
Discussion). Representative animals from China fell into nine 
clades that corresponded to nine species. Notably, A. uralensis 
from Xinjiang, China, fell into a clade (BS=100) comprised of A. 
pallipes from Xizang, China, and a sequence from GenBank 
(origin unknown), thus rendering A. pallipes paraphyletic 
(BS=69). A sole mitogenome representing “A. chejuensis" from 
Jeju Island was embedded in a clade containing A. agrarius. 
Apodemus draco, A. ilex, and A. semotus fell together in a 
well-supported clade (BS=100), but the relationships among 
the three species were not resolved (BS«50). Apodemus 
chevrieri, A. draco, A. ilex, A. latronum, and A. peninsulae also 
comprised subclades. 


Table 3 Measurement statistics of Rattus 


























SGL NBL ZB SBL 
n Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum 
H. nitidus 46 41.93 3.22 28.40 47.05 16.13 1.68 10.54 19.00 19.87 1.16 17.64 22.30 38.30 2.36 32.82 43.00 
R. nitidus nitidus 31 41.62 3.53 28.40 47.05 15.90 1.75 10.54 18.66 19.66 1.11 17.64 21.75 37.96 2.30 32.82 41.82 
H. nitidus thibetanus 15 42.56 2.44 38.39 45.85 16.62 1.47 14.02 19.00 20.31 1.15 18.63 22.30 38.99 2.40 35.19 43.00 
H. losea 31 37.12 3.32 31.24 46.53 13.51 1.77 10.06 18.35 17.99 1.38 15.02 21.12 34.37 2.61 29.72 40.55 
H. tanezumi 34 40.11 2.06 33.98 43.27 14.50 1.26 11.69 16.57 19.57 1.03 16.73 21.67 36.50 1.96 30.75 39.80 
H. andamanensis 44 41.49 3.06 34.87 46.80 15.13 1.75 11.25 17.84 20.33 1.63 16.88 24.05 37.69 2.93 31.76 42.72 
R. norvegicus 69 41.63 3.52 31.75 52.13 15.33 1.69 10.79 20.22 20.42 1.88 15.33 25.63 38.16 3.30 29.95 48.56 
R. exulans 4 32.86 2.15 31.09 35.88 11.33 1.64 9.91 13.14 15.77 0.73 15.15 16.79 30.49 1.22 29.15 32.04 
H. rattus 5 39.86 2.10 37.96 43.15 14.11 1.24 12.93 16.02 18.83 0.71 17.93 19.58 37.75 2.18 35.75 41.03 
UTRL UMRL ABL ML 
n Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum 
H. nitidus 46 21.29 1.38 18.62 24.36 7.39 0.54 6.15 8.92 7.36 0.56 6.03 8.32 21.70 2.02 16.41 26.17 
H. nitidus nitidus 31 20.96 1.25 18.62 23.02 7.29 0.53 6.15 8.35 7.28 0.59 6.03 8.32 21.08 1.85 16.41 24.65 
R. nitidus thibetanus 15 21.97 1.45 19.56 24.36 7.62 0.51 6.66 8.92 7.54 0.42 6.86 8.11 22.96 1.78 19.81 26.17 
H. losea 31 18.95 1.62 15.86 21.97 6.53 0.42 5.79 7.63 7.17 0.40 6.54 7.84 19.85 1.66 15.76 22.69 
H. tanezumi 34 20.13 1.22 16.76 22.47 7.38 0.38 6.39 8.18 7.45 0.51 6.22 8.62 20.76 1.09 17.81 22.59 
H. andamanensis 44 21.36 1.73 17.89 24.02 7.55 0.62 6.11 8.70 7.399 0.52 6.25 8.89 21.70 1.84 17.70 24.65 
R. norvegicus 69 21.66 1.93 17.01 27.14 7.35 0.47 6.32 8.55 7.40 0.56 6.51 8.68 22.02 2.33 17.02 29.63 
H. exulans 4 16.87 0.87 15.73 17.84 6.12 0.45 5.56 6.51 6.51 0.26 6.15 6.73 16.24 0.94 15.21 17.49 
H. rattus 5 19.61 0.83 18.83 20.88 6.73 0.38 6.41 7.39 7.41 0.48 6.99 8.24 19.48 0.98 18.43 20.70 
HBL TL HFL EL 
n Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum Mean SD Minimum Maximum 
H. nitidus 46 164.59 18.34 123.00 205.00 168.09 18.38 131.00 210.00 34.14 2.12 30.00 40.00 23.25 2.21 18.00 28.00 
H. nitidus nitidus 31 161.32 19.18 123.00 205.00 161.81 16.71 131.00 192.00 33.79 1.89 30.00 37.00 22.96 2.16 18.00 27.00 
R. nitidus thibetanus 15 171.33 14.84 155.00 205.00 181.07 14.81 162.00 210.00 34.87 2.45 31.00 40.00 23.87 2.26 19.00 28.00 
H. losea 31 152.07 20.82 103.00 192.00 145.58 15.64 110.00 192.00 28.58 2.04 24.00 32.00 19.12 2.95 12.00 27.00 
R. tanezumi 34 162.91 10.46 128.00 180.00 173.09 16.06 144.00 205.00 31.25 1.72 28.00 35.00 21.72 2.24 18.00 26.00 
R. andamanensis 44 173.02 19.54 126.00 212.00 186.89 24.45 115.00 231.00 31.65 1.87 27.00 36.00 23.15 2.29 19.00 32.00 
R. norvegicus 69 182.62 24.50 120.00 274.00 152.67 20.80 82.00 222.00 34.42 3.59 26.00 41.00 17.72 2.23 11.00 21.00 
H. exulans 4 115.50 1427 101.00 135.00 138.00 8.17 126.00 144.00 23.70 1.39 22.00 26.00 16.75 0.50 16.00 17.00 
H. rattus 5 170.80 16.99 150.00 194.00 190.60 11.80 180.00 210.00 33.00 2.24 30.00 36.00 22.60 1.82 20.00 25.00 


Abbreviations are explained in the Materials and Methods section. All measurements are in mm. 
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Figure 3 Principal component analysis of the first three principal components among the seven species of Rattus based on both 
external and craniomandibular variables (A) and craniomandibular variables only (B) 


Matrilineal genealogy of Rattus 


The interspecific relationships of Rattus using the mitogenome 
and cyt b sequences (n=396) were well-resolved (BS-95-100) 
or moderately resolved (BS-55-77) (Figure 5). Sequences 
representing animals from China fell into seven lineages that 
corresponded with A. nitidus, R. norvegicus, H. exulans, H. 
andamanensis, R. losea, R. rattus, and R. tanezumi. The clade 
of R. nitidus had two subclades, one from southern Xizang and 
the other from southeastern China (Figure 5). The tree depicted 
GenBank sequences deposited under different names within a 
shallow clade, most commonly with R. andamanensis. However, 
some specimens were also associated with R. losea, R. nitidus, 
R. tanezumi as well as R. nitidus from southern Xizang. 


DISCUSSION 


Genealogy and taxonomy 


Species of Apodemus are among the most destructive of 
all animal pests, yet little attention has been paid to their 
evolutionary relationships. Our trees were consistent with those 
from the robust study of Steppan & Schenk (2017), indicating 
the repeatability of both. However, the created molecular 
phylogenetic tree of Hattus was inconsistent with that of Aplin et al. 
(2011), which may be due to the different ways in which the trees 
were constructed (ML phylogeny here, but BI and NJ methods 
in Aplin et al. (2011)), different number of species, or different 
sequences of the cyt b gene (only two individuals of R. pyctoris 
(GenBank accession No. JN675511 and JN675512) from Aplin et 
al. (2011)). The unresolved relationships within Apodemus were 
not surprising and are likely due to early radiation in the evolution 
of this genus, as indicated by the saturation of the mitochondrial 
gene (Serizawa et al., 2000). Similar problems likely also occur in 
Rattus due to hybridization and introgression. Previously, Rattus 
was recovered as a paraphyletic genus (Steppan & Schenk, 201 7). 
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Fully resolved phylogenies require many slowly evolving and 
unlinked genes, which is not within the scope of the current study. 
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Figure 4 ML matrilineal genealogy of Apodemus derived from 
cyt b (Numbers above branches refer to bootstrap probabilities) 


Despite uncertainty in phylogenetic relationships, questions 
regarding taxonomy in both genera remain. The differences 
between A. pallipes and A. uralensis have been discussed 


previously in depth (Musser & Carleton, 2005). Our carefully 
identified specimens of A. pallipes were from southern Xizang 
(Pulan County). The average cyt b genetic distance between 
A. uralensis and A. pallipes was 5.4%, which was the smallest 
interspecific genetic distance in Apodemus. All our specimens 
of A. pallipes matched the original description and holotype 
(Musser & Carleton, 2005). Thus, A. pallipes undoubtedly 
occurs in China. The sequences of A. pallipes in GenBank 
were from Afghanistan and Pakistan, near the type locality of A. 
pallipes in Pamir Alta. However, as we had no access to these 
specimens, it was not possible to determine if they matched the 
morphological description of A. pallipes. 
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Johnson & Jones (1955) described A. chejuensis. Koh (1991) 
also recognized the species based on its large body size and 
mtDNA genotype. Corbet (1978) assigned it as a synonym of A. 
agrarius ningpoensis, whereas Musser & Carleton (2005) treated 
it as a synonym of A. agrarius. Our phylogeny embedded A. 
chejuensis in A. agrarius, and thus our results agree with the 
assignment of Musser & Carleton (2005). 

The taxonomic statues of A. draco remains uncertain. 
Apodemus ilex and A. semotus are close relatives to each 
other (Figure 4). Kaneko (2011) suggested that A. semotus did 
not differ significantly from A. draco. However, this endemic 
species of Taiwan was characterized by a dark gray pelage 
rather than the reddish-brown color of all other Asian species 
of Apodemus. Further, our ANOVA results demonstrated 
significant differences in TL and HFL between A. ilex and A. 
draco. Thus, we recognize all three as full species to better 


reflect their long evolutionary histories and distinct distribution 
patterns. Nevertheless, future comprehensive morphological 
diagnosis is desirable. 

Hodgson described R. pyctoris in 1845 from Nepal 
(Hodgson, 1845). This name was later replaced by R. rattoides 
or R. turkestanicus. Musser & Carleton (1993) resurrected the 
oldest name and it has been reported to occur in China (Allen, 
1940; Corbet, 1978; Ellerman & Morrison-Scott, 1951; Musser 
& Carleton, 1993, 2005; Wang, 2003). Feng et al. (1986) 
identified a series of specimens of R. pyctoris from Xizang and 
claimed that R. pyctoris closely resembled R. rattus but with 
a pale underbelly, relatively long nasal bone, and cusp t3 on 
M'. Our series of specimens from Xizang coincide with the 
characteristics of R. pyctoris described by Feng et al. (1986). 
However, phylogenetic analysis associated the species with H. 
nitidus. The original description and comments of Musser & 
Carleton (2005) on R. pyctoris point to its diagnostic characters 
as a very small cusp t3 on M', a wide and short rostrum 
(narrow and slender in R. nitidus), and chunky wide molars 
(thinner and gracile in H. nitidus). Except for the morphology 
of Mí, the Xizang specimens differed from R. pyctoris. 
Furthermore, many characters of the Xizang specimens 
also differed from R. nitidus, including the cusp t3 being 
present, gray-white underbelly, and larger measurements. The 
molecular phylogeny also placed the Xizang specimens and R. 
nitidus in different clades. Accordingly, we assign the Xizang 
specimens to a new, undescribed subspecies of H. nitidus. 

Peale (1848) described H. exulans from Society Island. 
Nevertheless, its existence in Taiwan, China has been 
recognized for a long time (Motokawa et al., 2001). The 
Guangdong Insects Institute collected specimens of R. exulans 
from Yongxing Island in 1975. Rattus exulans is the smallest 
Asian species in its genus. The specimens from Yongxing 
Island conformed to the characteristics of R. exulans. Thus, 
we confirm that R. exulans occurs in China in Yongxing Island 
and Taiwan. 

The earliest Chinese specimen of H. rattus (black type) was 
collected by A. B. Howell from Kuliang, Fukien in 1929 (Allen, 
1940). In 1955 and 1956, the Fujian Epidemic Prevention 
Station collected specimens from Fujian, which were confirmed 
by Shou (1962) as being R. rattus. Our examination of 
these specimens and one specimen from Guangdong Province 
resulted in the same conclusion. Thus, we confirm that H. 
rattus occurs in Fujian and Guangdong. 


Morphometrics- and molecular-based species identifications 


Regardless of skull and external measurements being similar 
between species, many interspecies measurements differed 
significantly. Species of Apodemus were easier to identify 
than Rattus. Furthermore, the different species of Apodemus 
exhibited stronger geographic distribution. | For example, 
although measurements could not discriminate between A. 
draco and A. ilex (current study) or A. semotus (Kaneko, 2011), 
all three were found to be allopatric: A. ilex occurs in Hengduan 
Mountains, south of the Yangtze River and west of the Jinsha 
River; A. semotus occurs in Taiwan only; and A. draco occurs 
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in the middle and lower reaches of the Yangtze River and in 
eastern China. Apodemus chevrieri, A. draco, and A. latronum 
co-occur in western Sichuan, but they were separated by the 


third upper molar and certain measurements (Figure 2A, B). 


Only one sequence of A. peninsulae in GenBank was likely 
misidentified (assuming no other error). Thus, the confusion 
between A. draco and A. ilex appears to be due to out-of-date 
taxonomy rather than misidentification. 

Identification of Rattus species using either morphometrics 
or molecular data requires caution. Unlike for the species 
of Apodemus, most species of Hattus are invasive in 
China and have likely experienced strong selection resulting 


in morphological modification to adapt to local habitats. 


Notwithstanding, it was possible to identify some species based 
on morphology alone, such as, R. andamanensis, which has a 
unique white belly, R. norvegicus, which has very short ears, 
and H. nitidus and R. norvegicus, which do not have the cusp 


t3 on M', with the former also having distinctly larger ears. 


The Chinese population of R. rattus is black all over its body, 
whereas H. exulans only occurs in islands of the South China 
Sea, including Taiwan, and has a very small head and body 
length. However, H. losea and H. tanezumi occur sympatrically 
in southern China. They are easily confused due to similar 
appearances and overlapping measurements. Most species 
showed significant overlap in the PCA plots (Figure 3A, B). 
Perhaps due to challenges in identification, GenBank contains 
many misidentifications. For example, sequences under the 
name of H. norvegicus occur in almost all clades (Figure 5). 

Our new sampling and survey of sequences supported the 
occurrence of nine species of Apodemus and seven species of 
Rattus in China. However, it is necessary to be cautious with 
morphometric and molecular analyses for species identification 
due to considerable intraspecific variation and considerable 
errors in GenBank. 


Alpha diversity of Apodemus and Rattus 

We determined that A. agrarius, A. chevrieri, A. draco, A. ilex, 
A. latronum, A. pallipes, A. peninsulae, A. semotus, and A. 
uralensis occur in China. In addition, considerable intraspecific 
diversity occurs in several species. Future comprehensive 
and integrative analyses can determine if further splitting is 
necessary and/or desirable. 

We determined that R. andamanensis, H. exulans, R. losea, 
R. nitidus, H. norvegicus, H. rattus, and R. tanezumi occur 
in China. Future research into the occurrence of H. pyctoris 
in China is not necessary. A new subspecies of H. nitidus is 
described as follows: 


Subspecies description 


Rattus nitidus thibetanus subsp. nov 

Holotype: Adult female, collected by Liao Rui on 15 January 
2011. The specimen was prepared as a skin with cleaned skull 
and deposited in the Sichuan Academy of Forestry (MT11197). 


Type locality: Motuo County, Xizang, China, N29.24344? and 
E95.169920°, 783 m a.s.l.. 
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Measurements of holotype: Weight: 179.6 g; HBL: 205 mm; 
TL: 200 mm; HFL: 40 mm; EL: 23 mm; SGL: 45.49 mm; SBL: 
43.00 mm; ZB: 19.98 mm; MB: 17.06 mm; ABL: 8.08 mm; 
LMxT: 7.38 mm; NBL: 19.00 mm. 


Paratypes: 5 specimens, with skins and skulls: XCY01001, 
°, 28.5048, 97.01045; XZ16259, ?; 27.47033, 88.91450; 
XZ16258, 3; 27.47033, 88.91450; MT020, 3; MTO035, d, 
29.25491, 95.21331. 


Additional specimens: 15 specimens (9 juveniles, 6 adults 
with skulls broken). Adults: XZ16260, °; XZ16253, 3; XZ11280, 
9: XZ11262, 3d; XZ11177, 3; and XZ11173, 9. Juveniles: 
XZ11207, $; MT11174, $; XZ11176, °; XZ11207, 2; XZ11208, 
$; XZ11232, $; XZ11175, 2; XZ11196, 3; and XZ11028, ©. 


Geographic distribution: The new subspecies is recorded 
from Yadong, Motuo, Nielamu, and Jilong counties, southern 
Xizang, China. 


Etymology: The name is derived from the type locality, 
southern Xizang (Tibet), China. 


Diagnosis: Cusp t3 present on M! in first transverse loop, 
but very small; head and body relatively large; tail length 
usually larger than head plus body length; belly gray-white; 
transition between darker dorsal and lighter ventral pelage 
abrupt; dorsum of feet white, not glossy. 


Description: Summer pelage from neck to hip uniform 
brown-black. Ventral hairs with gray-black base and gray-white 
tip, transition between darker dorsal and lighter ventral pelage 
relatively abrupt. Dorsal and ventral tail uniform brown-black; 
hairs on dorsal and venter of feet white, not glossy. 


Skull sturdy (Figure 6), in dorsal profile straight and brain 
case flattened; highest point of skull in middle of parietal bone. 
Nasal broad anteriorly narrowing posteriorly. Posterior margin 
of nasals irregular and protruding in front of maxilla. Posterior 
and anterior of frontal broad, middle narrower. Interparietal 
broad, anterior part triangle-shaped and posterior margin 
arc-shaped (Figure 6). Interorbital and temporal ridges present. 
Zygomatic arches medium in size, front part slightly broader. 
Auditory bullae moderately sized. Incisory foramen broad. 
Mandibles medium-sized (Figure 6). 

Upper incisors medium in size vertically downward and 
orange. Molars rooted; 1st upper molar with three transverse 
dental loops, first dental loop with 3 cusps, t3 present but small; 
2nd upper molar with three transverse dental loops, first dental 
loop only on lingual cusp; 3rd upper molar with three transverse 
dental loops, first dental loop only on lingual cusp, third loop 
only single semicircle and second loop rectangular; mandibular 
condyle and coronoid process large, but lower molar same as 
in other species of Rattus. 


Habitat: Specimens were collected from an abandoned 
farmland, along the footpath of a rice field where highland 
barley was grown, forest edge, shrubland, surrounding a house, 
and salvage station. 





10 mm 





10 mm 5mm 


Figure 6 Skull of new subspecies of Rattus nitidus 


Comparison with other subspecies: Compared with R. n. 
nitidus, t3 of the first dental loop present in Rattus nitidus 
thibetanus subsp. nov (vs. t3 absent or just vestigial; belly 
gray-white, and transition between darker dorsal and lighter 
ventral pelage relatively abrupt in Rattus nitidus thibetanus 
subsp. nov (vs. belly gray-white or yellow-gray, and transition 
vague in H. n. nitidus); dorsum of feet white, not glossy in 
Rattus nitidus thibetanus subsp. nov (vs. dorsum of feet white 
and shiny pearl in H. n. nitidus). The independent sample 
t-test demonstrated significant differences in UTRL, UMRL, ML, 
and TL between R. n. nitidus and R. n. thibetanus. The 
K2P distance for R. n. thibetanus and R. n. nitidus was 
0.019, smaller than the smallest interspecies distance known 
in Rattus. 
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A new genus of Asiatic short-tailed shrew (Soricidae, 
Eulipotyphla) based on molecular and morphological 
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ABSTRACT 


Blarinellini is a tribe of soricine shrews comprised of 
nine fossil genera and one extant genus. Blarinelline 
shrews were once widely distributed throughout 
Eurasia and North America, though only members 
of the Asiatic short-tailed shrew genus Blarinella 
currently persist (mostly in southwestern China and 
adjacent areas). Only three forms of Blarinella have 


been recognized as either species or subspecies. 


However, recent molecular studies indicated a 
strikingly deep divergence within the genus, implying 
the existence of a distinct genus-level lineage. We 
sequenced the complete mitochondrial genomes 
and one nuclear gene of three Asiatic short-tailed 
and two North American shrews and analyzed 
them morphometrically and morphologically. Our 
molecular analyses revealed that specimens ascribed 
to B. griselda formed two deeply diverged lineages, 
one a close relative to B. quadraticauda, whereas 
the other—comprised of topotype specimens from 
southern Gansu—diverged from other Blarinella in 
the middle Miocene (ca. 18.2 million years ago (Ma), 
95% confidence interval=13.4—23.6 Ma). Although 
the skulls were similarly shaped in both lineages, we 
observed several diagnostic characteristics, including 
the shape of the upper P^. |n consideration 
of the molecular and morphological evidence, we 
recognize B. griselda as the sole species of a new 
genus, namely, Pantherina gen. nov. Interestingly, 
some characteristics of Pantherina griselda are more 
similar to fossil genera, suggesting it represents 


an evolutionarily more primitive form than Blarinella. 


Science Press 


Recognition of this new genus sheds light on the 
systematics and evolutionary history of the tribe 
Blarinellini throughout Eurasia and North America. 


Keywords: Blarinellini; ^ Capture-hybridization; 
Mitogenome; Molecular phylogeny; Next-generation 
sequencing; Pantherina 


INTRODUCTION 


Asiatic short-tailed shrews, currently classified as species in the 
genus Blarinella, are small insectivorous mammals distributed 
mainly in central and southwestern China, adjacent Myanmar, 
and northern-most Vietnam. These small- to middle-sized 
shrews are uniformly black or dark brown and have large 
incisors, heavy tooth pigmentation, and a short tail that is 
typically 4096—6096 of the head-body length. The fore claws are 
enlarged, suggesting adaptation for a semi-fossorial lifestyle 
(Wilson & Mittermeier, 2018). 

The taxonomy of Blarinella has been studied since the 
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late 19'" century. The first recognized species was Sorex 
quadraticauda, described by Milne-Edwards & Milne-Edwards 
(1872) based on a specimen from Baoxing (=Mouping), 
northwestern Sichuan, China (the same type locality as 


the giant panda, Ailuropoda melanoleuca (David, 1869)). 


Milne-Edwards & Milne-Edwards (1872) documented the 
shrew's relatively short and somewhat square-shaped tail, 
well-developed incisors, and intensively dark pigmentation on 
the teeth. This species typically has five upper unicuspids, 
although the holotype specimen has only four (with the fifth 
one missing) on one side of its skull and five on the other, 
as discussed by Thomas (1911). When Thomas (1911) 
examined new specimens from Mt. Emei (Omi-San, 100 km 
south of Baoxing) in western Sichuan, S. quadraticauda was 
determined to be more closely related to the North American 
short-tailed shrews Blarina and Cryptotis (tribe Blarinini), than 
to the Old World Sorex, and thus assigned to its own genus 
Blarinella (literally “small Blarina"). 


Thomas himself recognized additional two Blarinella species. 


One (Blarinella griselda) was based on specimens from Lintan 
(=Taochou), Gansu, which were differentiated by their smaller 
size, grayer color, and shorter tails (Thomas, 1912). Thomas 
(1915) later recognized specimens from Pianma (=Hpimaw), 
Yunnan, as a third species, Blarinella wardi, based on their 
small size, dark color, and narrow skull. Since then, specimens 
collected across southern China and northern-most Vietnam 
have been assigned to one of these three groupings, though 
their taxonomic status has varied. For example, Allen (1938) 
and Hutterer (1993) placed B. griselda and B. wardi as 
subspecies in B. quadraticauda, whereas Corbet (1978) and 
Hutterer (2005) recognized all three as distinct species. 

Jiang et al. (2003) reviewed the taxonomy and distribution 
of the group and, based on multivariate and univariate 
morphometric analyses of skulls, recognized the three as 
distinct species but assigned only a few populations from 
northwestern Sichuan to B. quadraticauda. This three-species 
division has been widely accepted (Hutterer, 2005; Smith & 
Xie, 2008). Chen et al. (2012) was the first to apply molecular 
phylogenetic approaches to samples collected mostly from 
western Sichuan and Shaanxi. Their study revealed B. 
wardi as a distinct lineage at a basal position of the genus, 
confirming its species status, and found B. quadraticauda was 
a well-supported clade embedded within B. griselda, making 
the latter a paraphyletic group. The authors suggested that 
either this represented an incipient speciation event for B. 


quadraticauda, or the taxonomy of the genus was incorrect. 


It is worth noting that no holotype or topotype specimens 
of B. griselda from southern Gansu were included in either 
study, so the populations from Chongqing, Hubei, and eastern 
Yunnan were tentatively assigned to B. griselda based on their 
intermediate size between B. quadraticauda and B. wardi. 

All the studies mentioned above examined Blarinella 
taxonomy at the species/subspecies levels, implicitly accepting 
the monophyly of the genus. More recently, however, He 
(2011) discovered two genetically distinct, but sympatrically 
distributed, lineages from Mt. Qinling, Shaanxi, China, that 
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called this into question. One lineage clustered with previously 
sequenced Blarinella griselda, whereas the other formed a 
cryptic lineage that appears to have diverged from other 
Blarinella more than 10 million years ago (Ma), suggesting that 
an undescribed genus may exist. However, because the Mt. 
Qinling and southern Gansu mountains (including Lintan, the 
type locality of B. griselda) are on the same tectonic Qinling 
belt, and because no specimen from Gansu has ever been 
included in morphometric or genetic study, it was uncertain 
whether the cryptic lineage represented B. griselda or an 
unrecognized taxon. Despite efforts to explore the southern 
Gansu mountains (see He et al., 2013), no Blarinella species 
have been captured, and the question remains unresolved. 
Based on sequences from a single specimen from Gansu and 
a handful of specimens from Vietnam, Bannikova et al. (2017) 
confirmed the existence of two clades in the genus Blarinella. 
Because their sample locality in Gansu was near the type 
locality of B. griselda, their specimen is likely to be the “true” 
B. griselda. Based on the results of these two studies, we 
hypothesized that B. griselda occurs only in Gansu and Shaanxi. 
We further hypothesized that this distinct evolutionary lineage 
may represent a separate genus, and that ^B. griselda" from 
central and southern China, northern Vietnam, and Mt. Qinling 
is likely to be more closely related to B. quadraticauda. Thus, a 
revision of the taxonomic status of B. griselda is warranted. 

The systematic position of Blarinella has been revised 
by paleontologists based on craniodental characteristics and 
revisited by mammalogists using molecular data. Repenning 
(1967) firstly included the genus in the tribe Soricini, despite 
acknowledging the similarity between Blarinella and Blarinini 
(North American Blarina and Cryptotis). Reumer (1998) 
established the new tribe Blarinellini to include the extant 
Blarinella plus eight fossil genera. A close relationship between 
Blarinini and Blarinellini is supported by molecular phylogenetic 
analyses (Dubey et al., 2007; He et al., 2010). These lines of 
evidence together point to B/arinella as a relict genus of Blarinellini. 

The craniodental characteristics of Blarinella were first 
described by Repenning (1967) based on B. quaaraticauda 
and fossil species B. kormosi (-Zelceina kormosi; Storch, 
1995). Several characters recognized by Repenning (1967) 
were adopted by Reumer (1998) for defining the tribe 
Blarinellini, such as breadth of the interarticular area on the 
mandibular condyle, reduced posterior emargination on the 
upper molariform teeth, presence of entoconid crests on M, 
and Ms», and reduction of the talonid of M3. 

In the current study, we sequenced five complete 
mitochondrial genomes and seven partial ApoB genes of two 
Blarinella and five Blarina specimens and also compared their 
craniodental morphology. We used these data to update the 
taxonomy of Blarinella, focusing on the distinctiveness of B. 
griselda from Gansu and Shaanxi. 


MATERIALS AND METHODS 


Sampling and experiments 
All specimens and tissue samples were obtained following 
locally approved animal care procedures under the auspices 


of authorized collection permits. Blarinella specimens from 
China were collected by the Mammal Ecology and Evolution 
Group of the Kunming Institute of Zoology, Chinese Academy 
of Sciences, from Mts. Ailao, Gongshan, and Qinling (Figure 
1). Blarina brevicauda was collected near Bird Lake in 
Nopiming Provincial Park, Manitoba, Canada, while conducting 
research under Manitoba Conservation Permit WB12563. A 
tissue sample of Blarina hylophaga was obtained from the 
National Museum of Natural History, Smithsonian Institution, 


Washington DC, USA (USNM 568994) under transaction No. 


2073785 (Supplementary Table S1). DNA was extracted using 


Vietnam 





a DNA extraction kit (Qiagen DNeasy Blood & Tissue Kits, 
China) or the phenol/proteinase K/sodium dodecyl sulphate 
method (Sambrook et al., 2001). For each specimen, we 
sequenced the complete mitochondrial cyt b gene and a 
nuclear ApoB gene segment using Sanger sequencing. The 
primers were developed in previous studies (Dubey et al., 
2007; He et al., 2010) and are given in Supplementary Table 
S2. Each PCR product was sequenced using both the sense 
and anti-sense primers and were assembled using Lasergene 
SegMan v7 (DNASTAR). 


WB. griselda 

[73 B. quadraticauda 

BU B. wardi 

I B. cf. quadraticauda (B. griselda) 


———  National/province boudnary 


Figure 1 Sample localities of specimens used for molecular analyses 


Blarinella cf. quadraticauda (green dots) represents localities for specimens previously identified as B. griselda. Numbers refer to sample localities provided in 


the supplementary file. 


We selected one sample per species for Blarina brevicauda, 
Blarina hylophaga, and Blarinella wardi, and one for each of 
the two distinct lineages of B. griselda for sequencing of the 
complete mitochondrial genomes (mitogenome(s) hereafter) 
using next-generation sequencing (NGS). We used two 
approaches to obtain the mitogenomes, that is, long-range 
PCR and cross-species hybridization capture. Our full protocol 
has been described and published separately (Chen et al., 
2018). In short, we first amplified the mitogenomes using 
Phusion High-Fidelity DNA Polymerase (New England Biolabs, 
Canada) with two pairs of primers (Supplementary Table 
S2) designed within conserved regions of 12s rRNA, 16s 
rRNA, COX3, and ND1. The amplicons (~8 000 and 10 
000 bp in length) were purified using Serapure magnetic 
beads (Rohland & Reich, 2012) sheared to small fragments 
using NEBNext dsDNA Fragmentase (New England Biolabs, 
Canada) and ligated with barcode adapters (NEXTflex DNA 
Barcodes for lon Torrent, BIOO Scientific, USA) using a 
NEBNext Fast DNA Library Prep Set for lon Torrent (New 
England Biolabs, Canada) to construct DNA libraries. After 


ligation, the libraries were purified using Serapure magnetic 
beads and size-selected using a 296 E-gel on an E-Gel 
Electrophoresis System (Invitrogen, Canada). Libraries within 
the 450—500 bp size range were selected and re-amplified 
using a NEBNext High-Fidelity 2X PCR Master Mix (New 
England Biolabs, Canada). We then purified the libraries 
using Serapure magnetic beads and measured the DNA 
concentration of the purified libraries using Qubit Fluorometric 
Quantitation (Thermo Fisher Scientific, Canada). 

In cases where the sample could not be amplified 
successfully using the primers, we used an in-solution 
capture-hybridization approach to enrich mitochondrial DNA 
(Horn, 2012; Mason et al., 2011). Briefly, this approach 
uses mitochondrial probes and DNA libraries constructed using 
genomic DNA to capture mitochondrial-like sequences. To 
make the probes, we first amplified and purified two amplicons 
of the Blarinella griselda mitogenome. We then measured 
the DNA concentration of each amplicon using a Nanodrop 
2000 spectrophotometer (Thermo Fisher Scientific, Canada) 
and mixed the two amplicons to ensure the amount of DNA was 
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in proportion to their relative lengths. The mixed DNA was used 
to make mitogenome probes (hereafter termed baits) using a 
Biotin-Nick Translation mix kit (Roche, Germany) according to the 
manual. The baits were stored at —20 °C before hybridization. 
We constructed 450—500 bp-size libraries for each sample 
from previously extracted genomic DNA (each with a unique 
barcode as described above) using a NEBNext Fast DNA 
Library Prep Set for lon Torrent (New England Biolabs, 
Canada). The libraries were re-amplified and purified (see 
above). We mixed the baits and each DNA library at a ratio 
of approximately 1:10 and then incubated them for 24—48 h 
at 65 °C (Chen et al., 2018). The enriched libraries were 


reamplified and quantified using a Qubit as described above. 


Finally, we pooled the samples to ensure that each had a 
similar amount of DNA. We sequenced the samples using 
a v318 chip with the lon Torrent Personal Genome Machine 
(PGM). 


Molecular data processes and analyses 
After sequencing, we first binned the samples based on their 


sample-specific barcode adapters and converted the raw reads 
to FASTQ format using the Torrent Suite v4.0.2 (Thermo 
Fisher Scientific, Canada). We filtered out reads shorter than 


60 bp and trimmed off adapter sequences before assembly. 


We used MIRA4 for de novo assemblies (Chevreux et al., 
1999), during which we used a strategy (technology-iontor) 
specifically suitable for ion torrent data to better resolve 
the homopolymer insertions-deletions problem (Bragg et al., 
2013). We also used the Geneious iterative approach (up to 
100 iterations) to map the reads to a reference mitogenome 
of Blarinella quadraticauda from Baoxing, Sichuan, China 
(GenBank accession No. KJ131179.1). These assemblies 


were conducted using Geneious v8.1 (Kearse et al., 2012). 


We repeated both the MIRA4 and Geneious assemblies at 
least twice for each sample, aligned the assemblies for each 
sample using MUSCLE (Edgar, 2004), and generated a 5096 
consensus sequence in Geneious v8.1. We also compared 
the assembled mitogenome with the cyt b sequence of the 
same sample conducted using Sanger sequencing to ensure 
the mismatch between the two was «1 bp per 1 000 bp for 
quality control. Finally, we annotated the mitogenomes using 
the annotation function in Geneious v8.1. We removed the 
repeat region of the D-loop region assuming the short NGS 
reads were not long enough to correctly assemble. 

For phylogenetic analyses, we downloaded available 
mitogenomes of 13 soricine and one crocidurine (Crocidura 
attenuata) shrew from GenBank. We also included 16S rRNA 
and cyt b fragments of Blarinella, most of which were collected 
in one previous study (Chen et al., 2012), for a mitochondrial 
gene tree estimation. We aligned the mitogenomes with partial 
16S rRNA and cyt b using MUSCLE and removed all tRNAs, 
the ND6 gene (which is on the light chain), and the D-loop 
region from the alignment. The remaining 13 375-bp alignment, 
which included all coding genes on the heavy chain and two 
rRNA genes, was used for maximum-likelihood (ML) tree 
estimation using RAxML (Stamatakis, 2014) and implemented 
in CIPRES (Miller et al., 2010). We subdivided the alignment 
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by genes (15 partitions) and employed a GTR+G model for 
each gene. We conducted rapid bootstrapping and searched 
for the best-scoring ML tree, allowing the program to halt 
bootstrapping automatically under an extended majority rule 
criterion (autoMRE). We also subdivided the alignment by gene 
and codon positions (for coding genes) into 38 data blocks and 
conducted an additional analysis for comparison. 

We estimated an ApoB gene tree to ensure the mitogenomic 
tree was not strongly affected by incomplete lineage sorting or 
mitochondrial introgression. To accompany our newly collected 
sequences, we obtained sequences from previous studies 
(Dubey et al., 2007; He, 2011) and downloaded sequences 
of 10 soricid shrews as outgroups (Supplementary Table S1). 
We estimated the RAxML ApoB gene tree as described above. 
The ApoB gene was considered as a single partition because 
few mutations were observed in Blarinella sequences of the 
alignment. 

We estimated lineage divergence times from the 
mitogenome data. While we recognize that mitochondrial 
genes may overestimate true divergence times, this effect 
can be mitigated using appropriate calibrations. We 
applied a second calibration following Springer et al. (2018), 
which focused on the divergence times of eulipotyphlans, 
and estimated that the divergence between Crocidurinae 
and Soricinae occurred 36 Ma (95% confidence interval 
(Cl)=28.6—44.0 Ma). This is much older than fossil calibration 
used in previous studies (Dubey et al., 2007; He et al., 2010), 
which was based on the assumption that Crocidosoricinae 
was the ancestor of Crocidurinae, Myosoricinae, and Soricinae. 
Crocidosoricinae is now recognized as a tribe of Myosoricinae 
(Crocidosoricini), thereby invalidating that assumption (Furio 
et al., 2007). Fossils of both crocidurines and soricines are 
known from the Oligocene («34 Ma), with the oldest soricid 
fossil, Soricolestes soricavus, found from Middle Eocene strata, 
Khaychin Formation (Lopatin, 2002). Thus, the estimated 
divergence time in Springer et al. (2018) is congruent with fossil 
records. We estimated divergence times using BEAST v2.5 
(Bouckaert et al., 2014). We first estimated the best partition 
scheme and evolutionary model for each partition using 
PartitionFinder v2.1.1 (Lanfear et al., 2017). A four-partition 
scheme (Supplementary Table S3) was selected and the 
GTR+G model was supported as the best-fitting model for 
each partition. The monophyly of Soricinae was constrained 
so that Crocidura was fixed at the root of the tree. We employed 
the Birth-Death model as the tree prior and relaxed lognormal 
as the clock model prior. We used lognormal distribution for the 
prior of the divergence time of Soricidae. We set the mean to 
36 Ma, with a standard deviation of 0.135, so that the median 
of the prior was 35.7 Ma and the 95% CI was 28.6—44.5 Ma. 
We set 4x10? generations for each analysis and sampled 
every 4 000 generations. Convergence was assessed using 
Tracer v1.7 (Rambaut et al., 2018). We repeated the analyses 
twice and combined the tree files using LogCombiner v2.5 
and estimated the maximum clade credibility tree using Tree 
Annotator v2.5, both in the BEAST v2.5 package. 


Specimen examination and morphometric analyses 

We recorded morphological measurements from five newly 
collected Blarinella specimens from Shaanxi (representing the 
two distinct lineages of B. griselda), and included specimens 
used in our previous study, which were mainly deposited 
in the American Museum of Natural History, Smithsonian 
National Museum of Natural History, and Kunming Institute 
of Zoology (Appendix |; Jiang et al., 2003). External 
measurements, including head and body length (HB), tail 
length (TL), and hind foot length (HF), were recorded from 
specimen labels or field notes. Nineteen craniomandibular 
variables were measured using a digital caliper graduated to 
0.01 mm: condyloincisive length (CIL), interorbital breadth 
(IOB), cranial breadth (CB), cranial height (CH), maxillary 
breadth (MB), rostral length (RL), postrostral length (PRL), 
palatoincisive length (PIL), postpalatal length (PPL), upper 
toothrow length (UTR), maximum width across upper second 
molars (M2-M2), postpalatal depth (PPD), mandibular length 
(ML), lower toothrow length (LTR), length of lower incisor 
(LLI), zygomatic plate breadth (ZP), condyle-glenoid length 
(CGT), and breadth of coronoid process (BCP). We analyzed 
morphometric variation using craniomandibular variables for 
88 intact skulls by principal component analyses (PCA). We 
log4o-transformed each variable prior to conducting the PCA in 
SPSS v17.0 (SPSS Inc, Chicago, Illinois, USA). Based on the 
results of our molecular analyses (see Results), we assigned 
the single individual from Gansu and three of the Shaanxi 
specimens to B. griselda; the two other Shaanxi specimens 
previously identified as B. griselda were tentatively identified 
as B. cf. quaaraticauda. 


Morphological comparison and diagnosis 

We examined the morphology of the specimens as per 
Jiang et al. (2003), Repenning (1967), Reumer (1984), and 
Storch (1995), and followed their terminology for morphological 
descriptions. We took photos of skulls and teeth using a digital 
microscope VHX-2000C (KEYENCE). 


RESULTS 


Phylogenetic relationships 

Regardless of the partitioning scheme, each of our 
mitochondrial datasets revealed the same relationships with 
similar bootstrap values (BS); thus, only the tree partitioned 
by genes is presented (Figure 2). Nectogaline (Episoriculus, 
Neomys, and Nectogale) + anourosoricine (Anourosorex) 
shrews were strongly supported as a clade (BS=88) placed 
sister to the Soricini with moderate support (BS=75). Species 
of the genus Episoriculus were recovered as paraphyletic 
(BS=87). The Blarinellini (Blarinella) and Blarinini (Blarina) 
were strongly supported as a clade (BS-100), but, within 
Blarinellini, Blarinella formed two distinct sister clades (BS-95), 
the branch lengths of which were as long as those between 
Anourosoricini and Nectogalini. One of the two Blarinella 
clades contained the three specimens from Shaanxi and the 


single specimen from Gansu that we identified as B. griselda. 


Within the second major Blarinella clade, four specimens of B. 


wardi from northwestern Yunnan (Gongshan and Zhongdian) 
formed a distinct branch (BS=100) sister to the subclade 
consisting of the remaining specimens from Sichuan, Shaanxi, 
Yunnan, Chongqing, Hubei, and North Vietnam, previously 
identified as either B. quadraticauda or B. griselda (BS=100). 
There were several strongly supported subdivisions in the latter 
subclade. For example, samples from southern Shaanxi and 
northeastern-most Sichuan (Qingchuan, adjacent to Shaanxi; 
Figure 1) formed a basal position, sister to other members of 
this subclade (BS=73). Samples from Vietnam were placed in 
two distinct lineages (BS=100). 


In our ApoB nuclear gene tree, the same interspecific 
relationships were recovered in Blarinini-Blarinellini (Figure 3). 
A sister relationship between Blarinini and Blarinellini was 
strongly supported (BS=100), and two deeply diverged clades 
were recovered within Blarinella (BS=81). One comprised B. 
griselda from Shaanxi and Gansu. In the other clade, B. 
wardi again occupied a position (BS=97) sister to the remaining 
samples (i.e., B. quadraticauda). 


Our BEAST analyses recovered the same topology 
(Figure 4) as our mitochondrial and nuclear gene trees. 
All relationships were strongly supported with posterior 
probabilities >0.98, except for the sister-relationship between 
Episoriculus caudatus and E. macrurus, which was 0.90. 
The divergence time between Blarinini and Blarinellini was 
estimated to be 21.9 Ma (95% Cl=16.5—27.6 Ma), close to that 
estimated for Anourosoricini and Nectogalini (21.4 Ma, 95% 
Cl=19.2-30.8 Ma). The estimated divergence time between 
B. griselda and other Blarinella species was 18.2 Ma (9596 
CI213.4—23.6 Ma). The most recent common ancestor of B. 
quadraticauda and B. wardi was estimated at 5.8 Ma (95% 
Cl=3.8-8.2 Ma). 


Morphometric comparison 


External morphology and skull measurements are given 
in Table 1. Based on the PCA using craniomandibular 
measurements of intact skulls, the 15t principal component 
(PC1) accounted for 72.2% of the variation (eigenvalue=13.7), 
being positively correlated with all variables (loading >0.69), 
thus reflecting a size effect. The 2" principal component 
(PC2) accounted for 7.3% of the variation (eigenvalue=1.4) 
and was positively correlated with condyle-glenoid length 
(CGT), postpalatal length (PPL), and postrostral length (PRL), 
and negatively correlated with interorbital breadth (IOB), thus 
indicating strong correlation with the shape of the posterior 
cranium. On the PC1 and PC2 plot (Figure 5), B. griselda 
overlaps with B. cf. quaaraticauda specimens (i.e. those 
previously identified as B. griselda), suggesting they do not 
differ from each other by size or overall shape of the skull. 
When looking at the remaining specimens, B. quadraticauda 
occurs along the positive regions of PC1 and PC2, whereas 
B. wardi plots on the negative regions of both PC1 and 
PC2. Furthermore, B. cf. quadraticauda exhibits geographic 
variations. For example, the specimens from Chongqing 
are closest to B. cf. quadraticauda on the plot, being 
distinguishable from the specimens from Yunnan, Sichuan, and 
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Shaanxi. In addition, B. wardi from Yunnan and northern 
Myanmar also diverge, although they still overlap on the plot. 


Morphological diagnosis 


Because of the deep genetic divergence between B. griselda 
and the other Blarinella taxa, it is necessary to describe the 
morphology at higher taxonomic levels. Blarinella griselda 
exhibits entoconid crests on the lower M, and M (Figure 





6A), which are missing in Blarinini; its mandibular condyle 
also has a broad interarticular area (Figure 7A), which differs 
from Anourosoricini and Nectogalini, both of which have a 
narrow interarticular area; B. griselda has large but not fissident 
upper incisors (Figure 7B), which differs from Beremendiini with 
strongly fissident upper incisors. Finally, the M3 of B. griselda 
has a highly reduced talonid (Figure 7C), which differs from 
Soricini and Notiosoricini. 
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Figure 2 Mitochondrial gene tree 
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Branch lengths represent substitutions per site. Numbers above branches are bootstrap values supporting the relationship. Asterisks indicate bootstrap value is 


higher than 90. Key to type of gene data: 18MtG: mitogenome; 42 2 gene:16S rRNA and cyt b; 3 cyt b: only cyt b gene. 


The craniodental morphology of B. griselda closely matches 
descriptions given by Reumer (1998) for the tribe Blarinellini (see 
Systematic Biology section) with a single exception regarding the 
position of the entoconid and shape of the entoconid crest on the 
lower M, and Ms, which are variable in Blarinella (Figure 6), as 
discussed in a previous study (Storch, 1995). 

Though B. griselda does not differ from B. cf. quadraticauda 
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morphometrically (see Figure 5), there are several characters 
of the skull and teeth that differentiate B. griselda from B. 
quadraticauda, B. cf. quadraticauda, and B. wardi. For 
example, the tip of the coronoid of B. griselda is bent anteriorly 
relative to its position in the other species (Figure 8). The 
lambdoid crest (dorsal view) of B. griselda is triangle-shaped, 
whereas in the other species the crest is arcuate, extending 


more anteriorly (Figure 8). The lower P4 of B. griselda is 
small, whereas P, extends anteriorly in other species (Figure 
9A). On the lower molars of the other species, the paraconids 
extend anteriorly and the protoconids extend buccally, so that 
the trigonid basin area is large. Conversely, the paraconids 
and protoconids of B. griselda are not extended, and thus 
the trigonid basin is smaller (Figure 9A). In B. griselda, the 
entoconid and metaconid on M, and M, are close to each other 
but well separated, while the entoconid crests are present but 
distinctly lower than these two cusps (Figure 6A). In the other 
species, the entoconid and metaconid are either close to one 
another, connected by a high crest (B. quadraticauda; Figure 
6B, C), or the entoconid is inconspicuous (B. wardi; Figure 


6D), although there are some exceptions in B. quadraticauda. 


On unworn upper incisors, the apex of B. griselda extends 
more anteriorly, whereas the apex extends downward in the 
other species (Figure 8). The upper P^ of B. griselda has 
a triangular occlusal outline (Figure 9C), whereas P^ in the 
other species has a quadrangular occlusal outline (Figure 
9D). These nuanced but stable characteristics distinguish B. 
griselda from B. quadraticauda, B. cf. quadraticauda, and B. 
wardi. In consideration of the deep genetic divergence and 
morphological distinctions, we elevate B. griselda to a new 
genus, as described below. 


Systematic biology 
Family Soricidae G. Fischer, 1814 


Subfamily Soricinae G. Fischer, 1814 


Tribe Blarinellini Reumer, 1998 


The tribe was defined by Reumer (1998) basically based 
"horizontal ramus of mandible 


on dental morphology: 


short and high, making the lower dentition compressed 
anteroposteriorly and giving the lophs and lophids a 
compressed W-shaped appearance, mandibular condyle large, 
with a broad interarticular area; coronoid spicule of mandible 
well developed; teeth heavily pigmented; upper incisor 
protruding but not fissident; upper molariform teeth with 
a reduced posterior emargination, showing a tendency to 
develop a continuous endoloph; occlusal surface of M, nearly 
square; Ms with a reduced talonid". The “lower molars with 
the entoconid close to the metaconid so that the entoconid 
crest is short and high" are variable in B/arinella (Figure 6), as 
observed in a previous study (Storch, 1995). 
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Figure 4 Divergence times estimated using BEAST based on mitogenome data 


Branch lengths represent time. Number below each node denotes posterior probabilities (PP). Asterisks indicate PP is higher than 0.98. Numbers above 


branches refer to divergence time in millions of years. 
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Table 1 External and craniomandibular measurements, including mean values, standard deviations (top line), together with range 


and sample sizes (bottom line) of species included in the current study 


































































































Variable B. griselda B. quadraticauda B. cf. quadraticauda B. wardi 
HB 66.33+0.47 72.21+4.01 65.19+4.32 64.043-2.77 
66-67; 3 65-81; 19 52-73; 47 60-69; 23 
TL 36.67+1.7 45.63+4.13 36.37--2.38 37.57+3.14 
35-39; 3 40-60; 19 31-39; 46 32-43; 21 
HF 11+0 14.53+0.99 11.04+1.53 11.83+0.52 
11-1153 13-16; 19 8.5-13; 46 10.5-13; 23 
CIL 19.52+0.34 21.32+0.18 19.97+0.42 19.15+0.38 
19.14-19.96;4 21-21.69;19 19.13—20.93; 47 18.54-19.91; 23 
IOB 4.33+0.13 4.45+0.15 4.39+0.2 4.04+0.09 
4.2—4.52; 4 4.13—4.66; 19 4.06—4.83; 47 3.89-4.33; 23 
CB 9.08-40.26 9.5+0.1 9.12+0.22 8.33+0.22 
8.85-9.48; 4 9.34—9.7; 19 8.57-9.63; 46 7.84-8.72; 22 
CH 5.52+0.52 5.38+0.15 4.98+0.22 4.58+0.2 
4.62—5.89; 4 5.11—5.82; 19 4.43—5.15; 46 4.12-4.93; 23 
RL 6.51+0.15 7.68+0.25 7.25+0.35 6.644-0.2 
6.32-6.73; 4 7.19-8.05; 19 6.79-7.96; 47 6.25-7.11; 23 
PRL 11.49-+0.28 12.59+0.2 11.81+0.31 11.58+0.29 
11.2-11.96;4  12.28-12.94; 19 11.26-12.51; 46 11.06-12.11; 23 
MB 5.78+0.18 6.14+0.12 5.91+0.15 5.39+0.15 
5.57-6.05; 4 5.91—6.36; 19 5.62-6.26; 46 5.04—5.77; 23 
PIL 8.62+0.19 9.49--0.09 8.93+0.24 8.42+0.22 
8.42-8.9; 4 9.33-9.66; 19 8.26—9.45; 47 7.93-8.78; 23 
PPL 9.15--0.29 9.95+0.17 9.28--0.23 9.15+0.25 
8.76-9.5; 4 9.54—10.25; 19 8.7-9.76; 46 8.67—9.63; 23 
UTR 8.34+0.16 8.74+0.11 8.34+0.24 7.63+0.24 
8.16-8.59; 4 8.37—8.88; 19 7.79-8.64; 47 6.98-7.97; 23 
M?-M? 5.25+0.11 5.51+0.13 5.29+0.12 4.81+0.11 
5.15-5.44; 4 5.34—5.77; 19 5.08—5.58; 47 4.56-5.03; 23 
ML 10.27+0.21 10.7+0.15 10.06+0.32 9.52+0.29 
10.04-10.61;4 10.37-11;19 9.39-10.86; 46 8.91—9.95; 22 
LTR 7.83+0.16 8.6+0.1 8.1+0.2 7.5+0.16 
7.66-8.07; 4 8.41—8.77; 19 7.55-8.47; 46 7.08-7.74; 22 
LLI 4.49+0.07 5.21+0.14 4.73+0.2 4.39+0.18 
4.4—4.6; 4 4.84—5.49; 19 4.2—5.04; 46 4.13—4.79; 22 
CGT 8.69+0.17 9.27+0.14 8.74+0.18 8.61--0.2 
8.54—8.96; 4 9.04—9.55; 19 8.36-9.14; 46 8.13-8.96; 23 
BCP 1.19+0.01 1.34+0.05 1.2+0.09 1.11+0.09 
1.18-1.2; 4 1.24—1.45; 19 1.03-1.4; 46 0.98-1.34; 22 
MTL 4.73+0.15 5.23+0.08 5.032-0.11 4.68+0.11 
4.58—4.92; 4 5.12-5.41; 19 4.8—5.29; 47 4.37—4.87; 23 
PPL 3.32+0.08 3.64+0.08 3.47+0.1 3.21+0.08 
3.23-3.44; 4 3.48—3.78; 19 3.29—3.65; 47 3.07-3.34; 23 
ZP 2.03+0.09 2.28+0.09 2.14+0.1 1.97+0.11 
1.91-2.14;4 2.16—2.44; 19 1.86-2.38; 47 1.7-2.2; 23 


B. cf. quadraticauda represents specimens originally identified as B. griselda. Abbreviations of variable names are explained in the Materials 


and Methods. All measurements are in mm. 


Pantherina He Kai, gen. nov. 


Type species: Pantherina griselda (Thomas, 1912) 


Etymology: The genus name derives from the feminine Latin 
noun Panthera (“panther”) + the diminutive suffix -ina, hence, 
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‘little panther”. The name indicates animals in this genus are 
small but as black and aggressive as a panther. Griselda is a 


feminine given name. 


Suggested common name: Panther shrew; SJ}. 


Included species: Type species only. 


Diagnosis: External morphology and overall shape of 
skull very similar to Blarinella. Tip of coronoid process 
slightly bends anteriorly, differs from Blarinella (continuing 
vertically). |. Lambdoid crest angular in dorsal view (more 
rounded in Blarinella). Five upper unicuspids (four upper 
unicuspids in Petenyia). lą bicuspulate, with minute third 
posterior cuspule present (tricuspulate in Alloblarinella and 
Zelceina). Lower P4 small, not extending anteriorly (larger 
and extending anteriorly in B/arinella). Trigonid occlusal outline 
of lower molars trapezoidal-shaped (V-shaped or U-shaped in 
Blarinella). Apex of upper incisors extends anteriorly (directly 
downward in Blarinella). Upper P^ with triangular occlusal 
outline; protocone low, forming antero-lingual corner on ventral 
outline (trapezoidal-shaped in Blarinella). 


Description: Skull and mandible: Skull stout and robust. 


Rostrum short. Braincase dome-shaped but low. Occipital 


bone small, lambdoid crest triangle-shaped from dorsal view. 
Entoglenoid processes well developed. Horizontal ramus high. 


Coronoid process broad and spatulate-shaped, tip of coronoid 
bends anteriorly, anterior edge of coronoid concave. Coronoid 
spicule high and strongly pronounced, external temporal fossa 
shallow. Articular facets of condyle close to each other, with 
angle of intersection approximately 45°, upper facet narrow and 
cylinder-shaped, interarticular area broad. 

Lower teeth: Lower incisor (l4) obviously bicuspulate, with 
minute third posterior cuspule present. On buccal side, 
posterior end reaches level of anterior part of M;. On lingual 


side, slender groove obviously presents along lų toward tip. 


First lower unicuspid very small, larger part of tooth squeezed 


between |, and P4. Occlusal outline of P4 bulbous-shaped. 


Trigonids of M, and M, trapezoidal-shaped. On My, M, 
entoconids high and close to metaconid, with two cusps 
connected by low entoconid crest. M; with trigonid unreduced 
in size but with talonid reduced to single, low, nearly conical, 
slightly blade-like hypoconid, with no trace of entoconid crest. 

Upper teeth: Upper I! noticeably extends anteriorly. Apex 
of upper I! moderately long and broad, spatulate, not fissident; 
talon lower than first upper unicuspid. Five upper unicuspids 
(U) present. From occlusal view, U'—U® decreases in size 
of total area in order. U! twice as large as U?, and U4 
and US extremely small. U'—U* visible in lateral view of 
skull. P^ with no posterior emargination; protocone small, 
forming antero-lingual corner in outline of tooth; hypocone 
absent, hypoconal flange extends to approximately level of 
metacone; occlusal outline of P^ triangular in shape. M! and 
M? similar in shape: protocone low, forming antero-lingual 
corner in ventral outline, hypocone absent, hypoconal flange 
extends posteriorly; occlusal outline square-shaped. M? small, 
with well-developed paracone, lingual part consists of basin 
surrounded by U-shaped ridge. 


Remarks: Some characters observed in P griselda are more 


similar to fossil genera related to extant Blarinella species. 


For example, the anterior projection of the tip of the coronoid 
process in P griselda has also been observed in Petenyia 


hungarica and Zelceina kormosi (Reumer, 1984; Storch, 1995). 
The entoconid and entoconid crest on the lower M; and M, of P 
griselda have been observed in Petenyia hungarica (Reumer, 
1984). Of note, this is a variable character in Blarinella, with 
some B. quadraticauda also having a similar structure. The 
triangular-shaped P^ of P griselda is somewhat similar to that 
of Alloblarinella sinica and Petenyia hungarica and quite similar 
to that of Petenyia bubia (Reumer, 1984). 


Pantherina griselda 
Blarinella griselda (Thomas, 1912) 
Blarinella quadraticauda griselda (Allen, 1938) 


Type locality: 68 km SE Taochou (=Lingtan), 3 048 m a.s.l., 
Gansu, NW China. 


Holotype: Natural History Museum (British) No. 12.8.5.23, 


collected by J.A.C. Smith 

Suggested common name: Gray panther shrew; 17x 1f 
Measurements: Table 1 

Figure 6A, Figure 7, Figure 8 left, Figure 9A, C, E 


Distribution: Currently known from two localities in southern 
Gansu and southern Shaanxi (Figure 1). 


Description: As for the genus. 
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Figure 5 Plot of PC1 and PC2 scores from principal 
component analyses (PCA) of 19 log;9-transformed cranial 
measurements from 87 Blarinella specimens 

Most specimens were previously examined by Jiang et al. (2003). Specimens 
identified as B. cf. quadraticauda were originally identified as B. griselda in 
Jiang et al. (2003). 
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(A) B. griselda 


(C) B. quadraticauda 


0509390 


0905313 





0509420 


(B) B. cf. quadraticauda 





(D) B. wardi 


0611010 


Figure 6 Buccal view of lower M; and M, of four Blarinellini specimens 


Interspecific variations among entoconids and entoconid crests. The arrows point at the entoconids of M;. In B. wardi, entoconid does not present. 





Figure 7 Craniodental characters of B. griselda 


Articular facets of mandibular condyle (A), first upper incisor (B), lower Mg (C), buccal view (D) and lingual view (E) of lower incisor. 


Comparison: Specimens of Pantherina griselda from Shaanxi 
are superficially similar to Blarinella in gross morphology; 


however, a number of characters distinguish the two taxa. 


In P. griselda, the lambdoid crest is angular in dorsal view, 
whereas in species of Blarinella, the occipital bone extends 
anteriorly, resulting in a more rounded lambdoid crest (Figure 
8B). The tip of the coronoid process projects anteriorly rather 
than continuing vertically, as in Blarinella. The groove on the 
buccal side of |; is narrow, rather than broad in Blarinella. In 
the labial view, P4 of P griselda does not extend anteriorly, 
so the tooth and trigonid basin are smaller than that of 
Blarinella, in which the P4 extends anteriorly. In P griselda, 
the paraconids of M; and M» do not extend anteriorly, the 
trigonid occlusal outline is trapezoidal, and the trigonid basin 
area is small, whereas in Blarinella, the paraconids of M; and 
M; extend anteriorly, the metaconid extends posteriorly, the 
trigonid has a V-shaped or U-shaped occlusal outline, and the 
trigonid basin area is large. In P. griselda, the entoconid and 
metaconid on M; and M, are high and connected by distinctly 
lower entoconid crest.  Blarinella wardi has no entoconid 
or entoconid crest on M; or Mo, but the posterior faces of 
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its metaconids extend posteriorly forming a posterior ridge, 
which is not present in Pantherina. In B. quaaraticauda 
(including B. cf. quadraticauda), the entoconid crest is variable, 
depending on the locality. For example, some specimens 
from the type locality (i.e., northwestern Sichuan) have a 
low, moderately-developed entoconid located close to the 
metaconid, resulting in a short and moderately-high entoconid 
crest; in a few populations from northeastern Yunnan (B. cf. 
quadraticauda), the entoconid and entoconid crest are more 
similar to those of P griselda, although the height of the 
entoconid is lower; in other samples, the entoconid is not 
obviously present and is more similar to B. wardi. The upper 
I! of P. griselda extend anteriorly, with a long apex, whereas, 
in Blarinella, the apex is shorter and typically directed ventrally. 
The occlusal outline of P^ in P griselda is triangular and the 
protocone is small, forming an antero-lingual corner in the 
outline of the tooth, whereas in Blarinella, the occlusal outline 
of P^ is mostly trapezoidal, and the protocone is more robust, 
not forming an antero-lingual corner in the outline. On the M' 
and M? of P griselda, the hypoconal flange extends posteriorly, 
whereas in Blarinella, the hypoconal flange extends anteriorly. 





Figure 8 Mandibles and skulls of B. griselda (KIZ0509390; A,C,E,G) and B. quadraticauda (KIZ0905313; B,D,F,H) 


DISCUSSION 


Taxonomic implications 

Our results strongly support the conclusion that specimens 
ascribed to Blarinella griselda are paraphyletic; one is 
a close relative to B. quadraticauda, whereas the other 
presents an anciently diverged lineage of Blarinellini. The 
divergence time estimation (18.2 Ma) suggests a divergence 
more ancient than the ancestor of Nectogalini, and close to 
that between Anourosoricini and Nectogalini. Morphometric 
analyses could not differentiate the two forms, indicating 
that phenotypic evolution was not strongly associated with 
quantitative characters. However, several nuanced but stable 
diagnostic characters on the skull and dentition distinguish 
B. griselda from other species of Blarinella. Some of 
these craniodental differences were observed for several fossil 
genera of the tribe. In other words, the difference between 
B. griselda and members of the genus Blarinella is at the 


same level as that between currently recognized genera. 


Because the two extant lineages showed marked genetic 
divergence from one another (Figures 2, 3), and because of 
the observation of inter-generic level morphological variation 


(Figures 6—9), we recognize B. griselda as a distinct genus, 
namely Pantherina. 


It should be noted that Pantherina shares several characters 
with fossil taxa, such as a triangle-shaped P^ on occlusal 
view, indicating that the genus represents a more primitive 
form of the tribe and may have closer affinities with fossil taxa 
within Blarinellini. We have not attempted to fully revise the 
systematic paleontology in the current study, but taxonomic 
revision that includes fossils is warranted. We do not currently 
recognize Pantherina as one of the fossil genera based on the 
obvious differences in the skull and dentation. For example, 
it should not be assigned to Alloblarinella or Zelceina whose 
lower incisor is tricuspulate (bicuspulate in Pantherina). t 
should not be assigned to Petenyia which has only four 
unicuspids (Pantherina has five). Its small M3 also differs 
from that of Alloblarinella which has a large M; with complete 
trigonid and talonid. 


The new monotypic genus Pantherina is known currently 
only from southern Gansu (Lintan) and southern Shaanxi 
(Ningshaan). Despite extensive recent surveys (Chen et al., 
2012; K.H., unpublished data), members of this genus have 
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not been found in southern areas of China such as Chongqing, 
Hubei, or northwestern Sichuan—areas that are, however, 
occupied by B. quadraticauda. Pantherina is unlikely to be 
distributed in more northern areas due to absence of suitable 
habitat, and thus presumably has a restricted distribution in 
southern Shaanxi and southern Gansu. Because Pantherina 
has similar morphology to B. quadraticauda, and because they 
co-occur in Shaanxi, they might compete for the same habitat. 
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7 A 
- =~ a” ail 
Figure 9 Occlusal views of lower teeth and upper teeth, and 
occlusal view of upper P^ of B. griselda (A,C,E) and B. 


quadraticauda (B,D,F). Arrows point at the protocone of 
upper P^ 


We tentatively assign the remaining specimens previously 


identified as B. griselda from China and Vietnam to B. 


cf.  quadraticauda because they clustered together with 
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B. quadraticauda in a well-supported clade in both the 
mitochondrial (Figure 2) and nuclear gene trees (Figure 3), with 
the latter forming a lineage well embedded within this clade 
(Figure 2). This placement does not mean that a two-species 
scenario (B. quadraticauda and B. wardi) is sufficient to cover 
the species diversity of this genus. Instead, we believe that 
species diversity is still underestimated, as implied in our 
morphometric and phylogenetic analyses. For example, on the 
PCA plot, the animals from Chongqing are clearly distinguished 
from other B. quaaraticauda (Figure 5). The mitogenomic 
gene tree supports multiple fine-supported subclades within 
B. quaaraticauda (i.e., BS>70), which also show a strong 
geographic pattern. The animals from Shaanxi and 
northeastern-most Sichuan cluster together, forming a basal 
divergence within B. quadraticauda (BS=73). The specimens 
from Vietnam are also of two different originations (BS=100). 
Whether these distinguishable clades or morphometric groups 
represent undescribed species/subspecies or geographic 
populations remain open questions. 


Implications for systematics and macroevolution 


Recognition of a new genus of Blarinellini is especially 
interesting, not only because it is the second extant genus 
of the tribe, but because it could provide clues to fill 
the evolutionary relationship gaps among living and fossil 
taxa of this genera-rich tribe. Reumer (1998) assigned 
Blarinella together with eight fossil genera into Blarinellini, and 
Tregosorex, previously assigned to Blarinini, has also been 
assigned recently based on diagnosis of newly discovered 
fossil material (Doby, 2015). With the addition of Pantherina, 
11 genera can be recognized, which is the largest number 
for the Soricinae. The genus boundaries of Blarinellini 
shrews were established exclusively based on morphological 
diagnosis and comparison, and some species have been 
moved from one genus to another. The current classification 
and systematic hypotheses for the tribe Blarinellini are difficult 
to test using molecular-based approaches because there 
are no available resources. Although it is possible to 
use a morphological matrix, the results may be affected 
by homoplasy, especially when the number of characters is 
not large and molecular data are not integrated. Indeed, 
phylogenic analysis was carried out for Eurasian soricine 
shrews to include six genera of Blarinellini but revealed 
five polyphyletic lineages (Rofes & Cuenca-Bescos, 2009). 
Because all characters were equally weighted (which is usually 
not a good assumption for morphological characters), the 
topology was not constrained, and molecular data were not 
incorporated, high-level relationships among extant taxa were 
very different from the well-known molecular phylogenies: 
Blarinella was grouped with Anourosoricini (which is closely 
related with Nectogalini; see Figure 2), and Neomys (a 
nectogaline genus) was not grouped with the other Nectogalini 
shrews (see Figure 2). As acknowledged by the authors, the 
homoplasy index was high, so it may suffer convergence in 
morphology (Springer et al., 2013), which is very common 
in semi-fossorial small mammals (He et al., 2015). These 


previous studies suggest that a systematic relationship 
exclusively based on morphological data may be biased. The 
time-calibrated tree accomplished by a morphological matrix 
could provide a scaffold for analyzing craniodental evolution 
in the tribe, thereby assisting in examining hypotheses of 
systematic relationships. The latter could be especially true 
because Pantherina seems to represent a more primitive 
lineage than Blarinella, as some characteristics of Pantherina 
are also observed in fossil taxa. These include a more 
triangle-shaped P^, which is similar to that of Petenyia bubia 
from the Miocene. The presence of entoconids and entoconid 
crests on the lower M, and M, also seem to be pleiomorphic, 
which has been observed in many primitive fossil shrews, 


including Heterosoricinae from the Eocene (Repenning, 1967). 


Among the fossil genera of Blarinellini, four were exclusively 
distributed in Eurasia (Alloblarinella, Cokia, Hemisorex, and 
Paenepetenyia) and four were exclusively from North America 


(Alluvisorex, Anchiblarinella, Parydrosorex, and Tregosorex). 


In addition, Petenyia was distributed throughout both North 
America and Eurasia.  Finely-resolved phylogenetic and 
systematic relationships could help to illustrate the pattern of 
transcontinental migrations through time (Dubey et al., 2007). 
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Appendix I 
Specimens examined. 


Abbreviations: KIZ: Kunming Institute of Zoology, Chinese Academy of 
Sciences, AMNH: American Museum of Natural History, USNM: National 
Museum of Natural History, Smithsonian Institution. 


B. quadraticauda (n=19) 
Baoxin, Sichuan (KIZ820495, KIZ820588, KIZ820625, KIZ820626), 


Qionglai, Sichuan (AMNH111121, — AMNH111125, 
AMNH!1 11127, AMNH111128, AMNH111130, AMNH111136) 


AMNH111126, 


Wenchuan, Sichuan (AMNH111112, AMNH111114, AMNH111119, 
AMNH111137, AMNH111138, AMNH111140, AMNH111142, 
AMNH111145) 


B. cf. quadraticauda (identified as B. griselda in Jiang et al., 2003; n=45) 


Jingdong, Yunnan (KIZ640068, KIZ640105, KIZ640128, KIZ640129, 
KIZ640133, KIZ640134, KIZ640152, KIZ640191, KIZ640199, KIZ640200, 
KIZ640418, KIZ640429, KIZ640430, KIZ640431, KIZ640432, KIZ96180, 
KIZ98069, KIZ98107, KIZ98114, KIZ98151, KIZ98152, KIZ98313, 
KIZ98335, KIZ98378, KIZ98553, KIZ72053), Luchun, Yunnan (KIZ72053), 
Mt. Qinling, Shaanxi (KIZ0509420, KIZ0509423), Nanchuan, Sichuan 
(KIZ88126, KIZ88127, KIZ88252), Shimian, Sichuan (KIZ820742, 
KIZ820743, KIZ820744, KIZ820752, KIZ820753, KIZ820754, USNM574289, 
USNM574290, USNM574291, USNM574292, USNM574293), Wulong, 
Sichuan (KIZ88320, KIZ88321, KIZ88395) 


B. griselda (n=4) 


Gansu (AMNH60449), Mt. Qinling, Shaanxi (KIZ0509390, KIZ0509578, 
KIZ0509579) 


B. wardi (n=25) 


Bijang, Yunnan (KIZ820280), Degin, Yunnan (KIZ820280), Gongshan, 
Yunnan (KIZ73863, KIZ73889,  KIZ73910,  KIZ73920,  KIZ73927, 
KIZ73928, KIZ73929), Lijiang, Yunnan (USNM241403), Myanmar 
(AMNH1 14720, AMNH1 14721, AMNH1 14722, AMNH1 14723, 
AMNH1 14739, AMNH1 14740, AMNH114741, AMNH114748, 
AMNH1 14757, AMNH1 14758, AMNH114760), Pianma, Yunnan (KIZ74176), 
Weixi, Yunnan (KIZ810669), Yinjiang, Yunnan (KIZ76286), Zhongdian, 
Yunnan (KIZ810147) 


ZOOLOGICAL RESEARCH 


Taxonomic revision of the genus Mesechinus 
(Mammalia: Erinaceidae) with description of a new 


species 


Huai-Sen Ai'*, Kai He?**", Zhong-Zheng Chen?**, Jia-Qi Li”, Tao Wan”, Quan Li?, Wen-Hui Nie”, Jin-Huan Wang?, 


Wei-Ting Su”, Xue-Long Jiang?” 


! Gaoligongshan National Nature Reserve, Baoshan Yunnan 678000, China 

? Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming Yunnan 650223, China 

3 The Kyoto University Museum, Kyoto University, Kyoto 606-8501, Japan 

^ College of Life Sciences, Anhui Normal University, Wuhu Anhui 241000, China 

5 Nanjing Institute of Environmental Sciences, Ministry of Environmental Protection, Nanjing Jiangsu 210042, China 


ABSTRACT 


Hedgehogs in the genus Mesechinus (Family 
Erinaceidae), which include two currently recognized 
species (M. dauuricus and M. hughi), are distributed 
from northeast Mongolia to the upper Amur Basin in 
Russia and adjacent areas in northeast and northern 
China. In recent years, a population of Mesechinus 
hedgehogs was discovered from Mt. Gaoligong, 
southwestern Yunnan, China, far from the known 
distribution range of the genus. Furthermore, these 
hedgehogs are the only known population to be 
distributed at elevations higher than 2 100 m and in 
sympatry with gymnures. To evaluate the taxonomic 
status of these hedgehogs, we examined specimens 
representing Mesechinus taxa in China and further 


conducted morphometric and karyotypic analyses. 


Our results supported the existence of four species 
in China. Specifically, we identified the hedgehogs 
from Mt. Gaoligong as a new species, Mesechinus 
wangi sp. nov., and recognized M. miodon, previously 


considered as a synonym of either M. dauuricus or M. 


hughi, as a distinct species. Interestingly, we observed 
a supernumerary M^ on all specimens of Mesechinus 
wangi sp. nov., which is an extremely rare event in 
the evolution of mammalian dentition. 


Keywords: Mesechinus; Taxonomy; Morphometrics; 
Inhibitory cascade; Karyotype; New species; 
Supernumerary molar 
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INTRODUCTION 


Extant erinaceids, including spiny hedgehogs (Erinaceinae) and 
silky-skinned gymnures and moonrats (Galericinae), are found 
within the family Erinaceidae (Hutterer, 2005). The monophyly 
of each subfamily, as well as their sister-relationships, are 
well supported in various morphological and molecular studies 
(Corbet, 1988; Frost et al., 1991; He et al., 2012). These 
molecular studies also suggest that the living gymnures diverged 
from the ancestor of hedgehogs 40 million years ago (Ma), 
which is far older than the most recent common ancestor 
of living hedgehogs (Bannikova et al., 2014). Members 
in the two subfamilies are not only morphologically and 
genetically distinct but also characterized by different geographic 
distributions and habitats (Corbet, 1988). The living species in 
Galericinae are mainly distributed in humid montane forests of 
subtropical and tropical Southeast Asia (Echinosorex, Hylomys 
and Podogymnura), Southern China (Hylomys and Neotetracus), 


Received: 06 October 2017; Accepted: 02 April 2018; Online: 25 April 
2018 

Foundation items: This study was supported by the National Key 
Research and Development Program of China (2017YFC0505200), 
Fundamental Research Funds for NIES in 2017(GYZX170308) 
and Biodiversity Conservation Program by MEP of China, K. H. 
was supported by a JSPS Postdoctoral Fellowship for Overseas 
Researchers (P16092) 

* Authors contributed equally to this work 

"Corresponding authors, E-mail: hekai@mail.kiz.ac.cn; jiangxl mail. 
kiz.ac.cn 

DOI: 10.24272/j.issn.2095-8137.2018.034 


Zoological Research 39(5): 335-347, 2018 335 


and Hainan Island (Neohylomys). With their most recent 
common ancestor considered to be in the late Miocene 
(Bannikova et al., 2014), living hedgehogs have adapted 
to diverse habitats and are widely distributed throughout 
Africa (Atelerix and Paraechinus) and Eurasia (Erinaceus, 
Hemiechinus, Mesechinus and Paraechinus) in deciduous 
woodland, coniferous forest, forest steppe, grasslands, 
savanna, dry steppes, semi-desert, and even arid desert 
(Corbet & Hill, 1992); until recently, however, they have never 
been found in tropical or subtropical rainforest. In 2003, Ai 
(2007) discovered a small population of hedgehogs from the 
southern-most edge of Mt. Gaoligong in Yunnan Province at 
approximately 2 200-2 600 m a.s.l., near the border between 
China and Myanmar. These hedgehogs are characterized by 


the absence of a spineless section on their head and by ears 
of similar length to the surrounding spines, suggesting that 
they are members of the genus Mesechinus (Figure 1). The 
discovery was unexpected and of interest because: (1) the 
location is at least 1 000 km from the known distribution of any 
other hedgehog species; (2) the elevations are higher than that 
of any known hedgehog habitat; (3) the habitat is subtropical 
montane evergreen broad-leaved forest, which is typical habitat 
of the gymnures but differs from any known hedgehog habitat; 
and (4) the animals are sympatrically distributed with gymnures 
(Neotetracus sinensis), which is also the first ever record. While 
these clues indicate that the population represents a distinct 
taxon, its taxonomic status has yet to be resolved. 





Figure 1 Living Mesechinus wangi sp. nov. (KIZ 034115) 


Mesechinus hedgehogs are mainly distributed in northern 
China and Mongolia, as well as the Transbaikalia region 
and upper Amur Basin in Russia (Figure 2). Two species 
(M. dauuricus and M. hughi) were recognized in Mammal 
Species of the World (Hutterer, 2005). After Sundevall 
described the type species Erinaceus (Mesechinus) dauuricus 
in 1842, another five forms were recognized, including 
przewalskii Satunin, 1907, hughi Thomas, 1908, miodon 
Thomas, 1908, manchuricus Mori, 1927, and sylvaticus Ma, 
1964. Subsequently, however, manchuricus, przewalskii, 
and sibiricus were recognized as synonyms of M. dauuricus 
(Corbet, 1988), and miodon and sylvaticus as synonyms of 
M. hughi. The most debated species continues to be miodon, 
which was originally described together with hughi by Thomas 
(1908). Based on successive morphological research, some 
authors have included it in M. dauuricus (Corbet, 1978; Corbet 
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& Hill, 1992), whereas others have included it in M. hughi 
(Hoffmann & Lunde, 2008; Hutterer, 2005). Furthermore, 
karyotypic study of miodon from its type locality demonstrated 
variable chromosomal numbers ranging from 2n-44 to 48 
(Lin & Min, 1989; Kong et al., 2016a) due to the existence 
of B-chromosomes (Kong et al., 2016a), which has been 
interpreted as evidence of full species status (Kong et al., 
2016a, 2016b). However, B-chromosomes are rarely used 
for delimiting species and as its craniodental morphology has 
not yet been fully diagnosed, the species status of M. miodon 
remains suspicious. 

In this study, we integrated morphometric and karyotypic 
approaches to revisit the taxonomy of Mesechinus. We 
examined whether M. miodon is distinguishable from other 
species and were particularly interested in the taxonomic status 
of the hedgehog population found from Mt. Gaoligong. 


E 95° E 100° E 105° E 110° 


Russia 


Figure 2 Distribution of genus Mesechinus 


MATERIALS AND METHODS 


Specimens examined 

We examined 59 specimens (skins and skulls) of Mesechinus 
deposited in the Institute of Zoology (IOZ) and Kunming 
Institute of Zoology (KIZ) of the Chinese Academy of Sciences 
(CAS), Shaanxi Institute of Zoology (SXIZ), Northwest 
University (NWU), and China West Normal University (CWNU) 
(see Supplementary Appendix |). These specimens included 
the M. hughi sylvaticus holotype. Photo images of the M. 
miodon holotype were also obtained for examining diagnosable 
characters and for morphological description and comparison 
with other named species. Morphology was examined and 
described following Corbet (1988), Frost et al. (1991), 
Gould (1995), and Thomas (1908). Based on our diagnosis 
and comparison of external and craniodental morphology, 
we recognized four species/putative species, including M. 
dauuricus (to include M. dauuricus dauuricus (n=8) and M. 
dauuricus manchuricus (n-5)), M. hughi (to include M. hughi 
hughi (n=28) and M. hughi sylvaticus (n=3)), and M. miodon 
(n=9). We recognized the animals from Mt. Gaoligong as a 
new species, which we name herein as Mesechinus wangi sp. 
nov. (n=6). 
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Morphological measurement and analysis 


External measurements, including body weight (W), head and 
body length (HB), tail length (TL), length of hind foot (HF), and 
ear length (EL), were recorded from specimen tags. Spine 
length (SL) was measured from specimens. Twelve cranial 
characters were measured in millimeters with a digital caliper 
graduated to 0.01 mm (Table 1) following Pan et al. (2007): 
greatest length of skull (GLS), condylobasal length (CBL), 
basal length (BL), cranial height (CH), palatal length (PL), 
zygomatic breadth (ZMB), interorbital breadth (IOB), mastoid 
width (MTW), greatest width of nasal (GWN), breadth of first 
upper molar (BM'), length of upper tooth row (LUTR), and 
length of below tooth row (LBTR). We extracted measurements 
from Allen (1938) for the eight specimens of M. miodon 
deposited in the Natural History Museum. 

Morphometric variation was analyzed using principal 
component analysis (PCA) in SPSS v19.0 (SPSS Inc., Chicago, 
IL, USA). Only the 12 cranial measurements were used 
for PCA. All variables were log;o-transformed before PCA. 
One-way analysis of variance (ANOVA) was used to test 
significant differences in external and cranial variables among 
species. 
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Table 1 External and cranial measurements (mm) of Mesechinus specimens examined in this study (meantSD and range for each 
measurement and numbers of specimens measured (n) are given) 





Mesechinus wangi | Mesechinus miodon Mesechinus hughi Mesechinus dauuricus 
n=6 n=18* n=31 n=13 

W 411.20+48.66 505.00+168.73 341.39+127.82 562.41+130.37 
336.00-451.00;5  230.00-750.00; 6 112.00-750.00; 31 | 423.00-840.00; 11 

HB 202.40+26.10 195.22+24.26 189.71+24.20 206.21+22.30 
177.00—240.00; 5  120.00—220.00; 17  148.00-232.00; 31 175.00—261.00; 12 

TL 17.26+1.82 33.22+5.22 19.2343.32 24.08+3.65 
14.00-18.20; 5 25.00-43.00; 17 12.00-24.00; 27 17.00-30.30; 12 

HF 47.20+1.20 58.80+85.13 37.97+4.36 34.74+7.39 
45.30-48.00; 5 35.00—378.00; 16 30.00-47.00; 31 18.00-41.00; 12 

EL 29.60+1.74 28.81+3.13 22.94+3.99 31.19+3.44 
28.00-31.80; 5 24.00-34.50; 17 16.00-33.00; 31 22.30-34.00; 11 

GLS  54.75+0.81 54.10+2.18 49.39+1.58 55.18+3.21 
53.70—55.60; 4 49.30-57.20; 14 45.10-52.40; 23 50.20-58.40; 12 

CBL 54.55+0.68 53.18+2.47 48.46+1.61 54.72+2.94 
53.60-55.20; 4 48.50-56.30; 11 44.40-51.20; 23 49.40-57.40; 13 

CH 17.13+0.69 18.67+0.72 16.14+0.97 18.37+0.58 
16.10-17.60; 4 17.80-19.70; 6 14.90-18.20; 21 17.20-19.10; 9 

BL 50.00+1.58 49.64+2.12 45.55+1.32 51.83+2.02 
47.70—51.30; 4 44.70—52.30; 14 43.20-48.80; 21 48.10—54.50; 13 

PL 30.25+0.58 28.82+1.46 26.58+0.63 28.60: 1 
29.50-30.80; 4 27.00-32.18; 14 25.70-28.40; 21 

ZMB  33.97+0.23 32.7742.17 28.90+1.72 32.62+2.93 
33.70-34.10; 3 28.70-37.08; 14 25.70-32.00; 22 28.40-36.40; 13 

IOB 14.68+0.38 13.87+0.83 12.51+0.52 13.86+0.72 
14.20-15.10; 4 12.90-15.10; 6 11.70-13.60; 23 13.00-15.10; 9 

MTW  25.60+0.73 25.93+1.23 21.67+1.60 25.58; 1 
24.70—26.20; 4 24.30-28.30; 14 19.50-24.50; 21 

GWN  4.30+0.00 2.70+0.23 2.97+0.30 2.96: 1 
4.30—4.30; 3 2.37—2.94; 6 2.60—3.60; 23 

BM! 21.43+0.31 21.08+0.69 17.38+0.77 T 
21.10-21.70; 3 20.30-22.30; 14 16.50-19.50; 21 | 

LUTR 27.90+1.18 27.25+1.03 24.65+1.15 27.85: 13 
26.70-29.10; 4 25.70—29.02; 14 21.40—26.10; 23 

LBTR 24.85+0.51 24.91+0.73 21.19+0.80 24.30; 1 
24.20-25.30; 4 23.40-25.70; 14 20.20-23.70; 21 


*: Includes measurement of nine specimens measured by Allen (1938). Abbreviations are given in the Materials and Methods section. 


Cell culture and cytological preparation 


One specimen representing Mesechinus wangi sp. nov. 
(museum catalog number: KIZ 034115) was used for cell 
cultures. We followed Hungerford (1965) for cell culture and 
metaphase preparation. The fibroblast cell cultures derived 
from skin fibroblasts and bone marrow are stored in the 
Kunming Cell Bank, Kunming, Yunnan, China. Images were 
captured using the Genus System (Applied Imaging Corp., 
USA) with a CCD camera mounted on a Zeiss Axioplan 2 
microscope. Chromosomes of Mesechinus wangi sp. nov. 
were arranged based on their relative length in order from 
longest to shortest. 
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RESULTS 


Morphological comparison and diagnosis 

As mentioned previously, the hedgehogs from Mt. Gaoligong 
could be assigned to Mesechinus unambiguously based on 
external morphology. These animals lack a spineless area 
on their heads, which is distinct from Atelerix, Erinaceus, and 
Paraechinus (Figures 2, 3), and their ears are similar to the 
surrounding spines in length, which is distinguishable from 
Hemiechinus. The skull and teeth are also characterized by 
several typical Mesechinus features, including a robust jugal 
reaching the lacrimal (Figure 4), shallow suprameatal fossa, 
and narrowly separated anterior and posterior borders of the 


suprameatal fossa (Figure 4), which distinguish it from all 
other genera (see Frost et al., 1991 for discussion). We 
compared the external and craniodental morphology of our 
specimens. It is worth noting that the sample size for some 
species/subspecies was small and may not reflect intra-specific 
variation, especially that of teeth (see discussion in Frost et al., 
1991; Gould, 2001), which needs to be verified in future study. 

All specimens examined in the current study showed few 
wholly white spines (Corbet, 1988; Figure 3). Spine lengths 
from longest to shortest were: M. miodon (~26 mm), M. 
dauuricus and Mesechinus wangi sp. nov. (~21—24 mm), and 
M. hughi (^21 mm). Spine color pattern was used by Thomas 
(1908) as a distinguishing feature for describing M. hughi and 





Figure 3 External morphs and spines of Mesechinus wangi sp. 
hughi (KIZ 027029) (C), and M. d. dauuricus (KIZ 027005) (D) 


Based on skull morphology, frontals were relatively higher 
than parietals in M. hughi and Mesechinus wangi sp. nov., 
whereas parietals were higher than frontals in M. dauuricus and 
M. miodon (Figure 4). On the ventral side of the skulls, a 
posterior palatal shelf and well-developed spine were present on 
all specimens examined (Figure 4). Mesechinus miodon showed 
longer spines than other species (Figure 4B). An epipterygoid 
process was present in all specimens (Figure 4; Frost et al., 
1991). In M. miodon this process was well developed, extending 
labially (Figure 4B), but only slightly or moderately developed in 


M. miodon (Figure 3). We found that M. hughi from Shaanxi 
(including topotype of M. hughi hughi) and Shanxi (holotype of 
M. h. sylvaticus) and M. dauuricus shared similar characters: 
that is, white for two-thirds of length, followed by black ring, 
narrow light ring, and black tip (Figure 3C, D). Mesechinus 
miodon was distinguished by spine light brown (rather than 
wholly white) for two-thirds of its length, followed by broad 
blackish-brown rings (rather than wholly black), light brown 
terminal (3-4 mm), and non-black tip (Figure 3D). Mesechinus 
wangi sp. nov. was differentiated from M. hughi by dark ring 
extending to tip on most spines, with narrow white ring near tip 
(Figure 3A). 





nov. (Type KIZ 022028) (A), M. miodon (Type BM 9.1.1.9) (B), M. h. 


other taxa (Figure 4). The basisphenoid of M. dauuricus was 
previously considered to be uninflated, intermediate between the 
condition of Hemiechinus and that of Atelerix and Erinaceus 
(Frost et al., 1991). According to our examination, however, 
the basisphenoid was inflated in M. dauuricus, M. hughi, and 
M. miodon (Figure 4B, C), similar to the condition observed in 
Hemiechinus auritus, whereas the basisphenoid was uninflated 
in Mesechinus wangi sp. nov. (Figure 4A), similar to that of 
Atelerix and Erinaceus. 

On the dorsal side of the skull, the nasal-maxilla relationship 
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was used in Corbet (1988) and Ma (1964), though Frost 
et al. (1991) determined that the relationship exhibited too 
much inter-specific variation. | Nevertheless, nasal breadth 
was obviously and significantly (see below) different between 
Mesechinus wangi sp. nov. and other species. More specifically, 
Mesechinus wangi sp. nov. was characterized by: nasal 
broad, premaxilla extending only slightly posteriorly and frontal 
extending only slightly anteriorly on dorsal side, premaxilla not 
touching frontal, and nasal and maxilla sharing long common 
sutures (Figure 4A). All other species exhibited much narrower 
nasal (Figure 4B-D). Mesechinus hughi could be characterized 
by: premaxilla extending posteriorly, frontal extending anteriorly, 
not touching premaxilla, with nasal and maxilla sharing short 
sutures (Figure 4C). Mesechinus d. manchuricus and M. miodon 
could be diagnosed by: premaxilla extending posteriorly, frontal 
extending anteriorly, premaxilla and frontal touching on dorsal 
side of skull (or nearly so), with nasal and maxilla not sharing 
common suture (Figure 4B, D). 

As M. miodon was named based on its small P? (triangular 
(equal-sided) in shape; Thomas, 1908)) examination of 
teeth was unavoidable here. Mesechinus dauuricus could be 
diagnosed by P? similar to P? in size; M. miodon could be 
diagnosed by P? smaller than P? (and smaller than that of M. 
dauuricus); Mesechinus wangi sp. nov. and M. hughi could be 
diagnosed by P? small, though similar to M. miodon (Figure 4). 
All species showed reduced upper M? and small trigonid (Figure 
4). Most notably, Mesechinus wangi sp. nov. could be further 
distinguished by consistent presence of single-rooted M^ on all 
specimens examined (Figure 5), much smaller than M$. 


Morphometric analyses 
External and cranial measurements of each species are given 


in Table 1. Thirty intact skulls were used for PCA, including 
specimens of M. dauuricus manchuricus (n=1), M. hughi hughi 
(n=20), M. miodon (n=6), and Mesechinus wangi sp. nov. 
(nz3). 

The first principal component (PC1) accounted for 74.24% of 
variation (eigenvalue=8.91) and was positively correlated with 
all variables, reflecting a size effect (Table 2). The second 
principal component (PC2) accounted for 9.4096 of variation 
(eigenvalue=1.13) and was dominated by MTW (loading=0.93), 
but was also positively correlated with PL, BM', LBTR, BL, 
and LUTR (loading>0.53). The third principal component 
(PC3) represented 4.90% of variation (eigenvalue=0.59) and was 
correlated primarily with GWN (loading=0.97). 

As shown in Figure 6A, M. dauuricus, M. miodon, and 
Mesechinus wangi sp. nov. plotted closely in the positive 
region of PC1 and PC2, indicating that these taxa had larger 
skulls. Further, M. hughi plotted in the negative region of PC1, 
indicating this species had a smaller skull. In the PC1 and PC3 
figure (Figure 6B), Mesechinus wangi sp. nov. plotted in the 
positive region of PC3 against all other species, indicating this 
species had the widest nasal. 

We employed one-way ANOVA for all external and cranial 
variables. The results showed that all variables were 


significantly different among the four species (P<0.001), except 
for HB (F=2.080, P=0.134), HF (F=0.522, P=0.596), and PL 
(F=7.561, P=0.002). 





Figure 4 Dorsal ventral and lateral views of skull and mandible of Mesechinus wangi sp. nov. (Type KIZ 022028) (A), M. miodon 
(Type BM 9.1.1.9) (B), M. h. hughi (KIZ 027029) (C), and M. d. dauuricus (KIZ 027005) (D) 


340 www.zoores.ac.cn 





Figure 5 Right upper molars of Mesechinus wangi sp. nov. (Type KIZ 022028) (A), M. miodon (Type BM 9.1.1.9) (B), M. h. hughi 


(KIZ 027029) (C), and M. d. dauuricus (KIZ 027005) (D) 


Table 2 Factor loading eigenvalues and percentage of variance 
explained for PC1, PC2, and PC3 from principal component 








analysis 
; Component 
Variables 
1 2 3 

CH 0.916 0.200 —0.136 
CBL 0.836 0.461 0.221 
GLS 0.830 0.483 0.160 
BL 0.800 0.541 0.141 
LUTR 0.720 0.5388 0.236 
ZMB 0.632 0.450 0.288 
IOB 0.542 0.480 0.458 
MTW 0.293 0.928 0.026 
PL 0.523 0.721 0.292 
BM! 0.623 0.696 0.216 
LBTR 0.631 0.683 0.155 
GWN 0.034 0.093 0.970 
Eigenvalues 8.909 1.128 0.588 


Total variance explained (%) 74.241 9.403 4.899 
Abbreviations are given in the Materials and Methods section. 
Karyotypic characteristics of Mesechinus wangi sp. nov. 


The karyotypes of Mesechinus wangi sp. nov. are shown in 
Figure 7. The diploid number (2n) and autosomal fundamental 


number (FNa) were 48 and 92, respectively (Figure 7A). The 
autosomes and X chromosomes were biarmed; however, we 
could not determine whether the Y chromosome was biarmed 
as it was too small. In total, 22 metacentric + 24 submetacentric 
autosomes were found in the karyotype. Both the X and Y 
chromosomes were metacentric, with the Y chromosome being 
smallest. G-banded karyotypic analysis identified homologous 
chromosomes (Figure 7B). 

Compared with other species in the genus Mesechinus, 
Mesechinus wangi sp. nov. had the same 2n and FNa as 
M. dauuricus and M. hughi, but differed from the reported 
karyotype of M. miodon, which is characterized by the presence 
of 0—4 B-chromosomes (2n244—48; FNa-82-92; Table 3). The 
numbers of metacentric chromosomes (M), submetacentric 
chromosomes (SM), and subtelocentric chromosomes (ST) 
also differed among species. 


DISCUSSION 


Taxonomic implications 


We compared the morphology and karyotypes among 
Mesechinus taxa in China. Although sample sizes were small 
for several forms, the patterns detected in the morphological 
and karyotypic analyses helped clarify the taxonomy of this 
genus. lt is worth noting, as well discussed in Gould (2001), 
that dental structures in hedgehogs can exhibit considerable 
intraspecific variation, and all dental characters should be 
treated with caution. 
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Figure 6 Plot of Mesechinus spp. for principal component factors 1 and 2 (A) and 1 and 3 (B) 


Mesechinus miodon is still recognized as a subspecies 
of either M. dauuricus (Corbet & Hill, 1992) or M. hughi 
(Hutterer, 2005). Here, however, we recognized M. miodon 
as a distinct species based on morphometric and karyotypic 
analyses (Figure 6; Table 3). Thomas (1908) named M. 
miodon and M. hughi but did not compare either with M. 
dauuricus. |n the current study, M. miodon was easily 
distinguished from M. hughi (Figure 6A), with its obviously 
larger cranial measurements (Table 1). Mesechinus miodon 
and M. dauuricus exhibited similar overall skull shape and size 
(Table 1, Figure 6), but M. miodon was distinguishable based 
on different spine color pattern and smaller P?. 

The implications of the karyotypic evidence are two-fold. 
On the one hand, M. miodon had a smallest number of 
metacentric chromosomes in the genus, and the numbers of 
submetacentric and subtelocentric chromosomes were also 
different from that of M. dauuricus, indicating that these two 
morphologically similar forms were distinct species. On the 
other hand, the existence of B chromosomes in the topotype 
of M. miodon (from Yulin), as reported in previous studies (Lin 
& Min, 1989; Kong et al., 20162), should be treated with caution. 
Although B chromosomes are heterochromatic and can be 
verified easily using the C-banding karyotypic approach (e.g., 
Badenhorst et al., 2009), this was not adopted in the original 
studies mentioned above. The number of B-chromosomes is 
usually stable, rather than highly variable as reported for M. 
miodon (0—4), and is usually an odd number, rather than an 
even number (Table 3) as reported in other mammals (e.g., 
Badenhorst et al., 2009). Thus, reexamination of the C-banding 
karyotype using additional samples is warranted. Finally, 
B-chromosomes are considered adaptive characters that can 
vary between populations and may be a poor characteristic for 
distinguishing species. 
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Figure 7 Karyotypes of Mesechinus sp. (KIZ 034115) 


A: Conventional karyotype; B: G-banding karyotype. 


We recognized the hedgehogs from Mt. Gaoligong as a 
distinct new species due to their many unique features. This 
new species is the first known hedgehog to be found at 
elevations higher than 2 100 m (Erinaceus europaeus and 
M. hughi are distributed no higher than 2 100 m), while also 
inhabiting subtropical evergreen broad-leaved forests (Figure 8) 
and co-occurring with gymnures. The color pattern of its 
spine is distinguishable from other Mesechinus species due 


Table 3 Karyotypes of the four recognized Mesechinus species 


to the lack of a narrow white ring (Figure 3A). It has a broad 
nasal that shares a long common suture with the maxilla 
(Figure 4A), which differs from all other taxa. The presence of 
a supernumerary M* is also highly distinctive (Figure 5A). We 
propose the name of Mesechinus wangi sp. nov. for the new 
species, in memory of the late Prof. Ying-Xiang Wang, a highly 
respected mammologist from the Kunming Institute of Zoology, 
Chinese Academy of Sciences. 





Species 2n NFa'  Autosomes X/Y chromosomes Locality Reference 

M. dauuricus 48 92 22M+14SM+10ST SM/T Chita, Russia Korablev et al., 1996 

M. hughi 48 92 30M+12SM+4ST M/T Yulin, China Lin & Min, 1989 

M. miodon 44 82 18M+24SM ST/T Yulin, China Lin & Min, 19892 
44-48° 84-92 18M+24SM+0-4B(M/SM) ST/T Yulin, China Kong et al., 2016 

M. wangi sp. nov. 48 92 22M+24SM M/M4 Baoshan, China This study 


1: Diploid chromosomes classified into metacentric (M), submetacentric (SM), subtelocentric (ST), and telocentric chromosomes (T). ?: 


Animals were identified as M. dauuricus in the original article. Kong et al. (2016a) argued that the specimens should be M. miodon. ?: 
Because of the existence of 0-4 B-chromosomes (M or SM; Kong et al., 2016a), the 2n could be 44—48 and NFa could be 84-92. ^: Y 
chromosome most likely biarmed (see Figure 6B), but was too small to be confirmed. 





Figure 8 Habitats of Mesechinus wangi sp. nov. within Gaoligongshan National Nature Reserve 


Supernumerary upper molar 

During the evolution of mammalian dentition, changes in 
tooth number are very common between taxa and within 
species, including in the eulipotyphlan mammals. Differences 
in number of teeth can be used as diagnostic characteristics 
in shrews (e.g., between Crocidura and Suncus) and talpid 
moles (e.g., between Euroscaptor Mogera, and Parascaptor), 
especially the number of incisors and premolars (or unicuspids 
in Soricidae). These differences are often triggered by 
geographic isolation and speciation, such as observed in 
the Persian mole (Talpa davidiana; KryStufek et al., 2001) 
and Japanese mole (Mogera wogura; Asahara et al., 2012). 


Nevertheless, an increased number of molars is an extremely 
rare event observed in only a few taxa, such as the bat-eared 
fox (Otocyon megalotis), and can be impacted by the evolution 
of different feeding behaviors and explained using an inhibitory 
cascade model (Asahara, 2013; Asahara, 2016; Asahara et 
al., 2016). For example, adaptation toward an increased bite 
force in the hypercarnivorous bush dog (Speothos venaticus) 
resulted in an enlarged M4, thus prohibiting the development of 
M; dentition (Damasceno et al., 2013). Although erinaceids are 
characterized by highly variable dentition and tooth structure 
(Gould, 2001), M4 and M^ have never been observed in 
either living or fossil species. Loss of M; has been observed 


Zoological Research 39(5): 335-347, 2018 343 


in short-faced hedgehogs, an extinct subfamily of gymnures 
(Brachyericinae), following remarkable shortening of the skull 
(Rich & Rich, 1971). Based on teeth features, Lopatin 
(2006) hypothesized that brachyericines from Asia underwent 
adaptation toward carnivory. If the inhibitory cascade model 
is valid in Erinaceidae, the presence of a supernumerary M4 
in Mesechinus wangi sp. nov. could be attributed to a 
combination of genetic bottleneck and isolation (hypothetically 
after long-distance dispersal from northern China) as well 
as adaptive selection from an omnivorous-insectivorous diet 
toward a highly omnivorous one resulting in reduced inhibition 
during the development of the upper molars (Asahara, 2016; 
Asahara et al., 2016). 


Taxonomic accounts 
Mesechinus dauuricus (Sundevall, 1842) Daurian Hedgehog 


Erinaceus dauuricus Sundevall, 1842: 237. Type locality: 
Dauuria, Transbaikalia, USSR. 
Hemiechinus przewalskii Satunin, 1907: 181. Type locality: 


North China. 

Hemiechinus manchuricus Mori, 1927: 108—109. Type locality: 
"Koshurei, South Manchuria” (=Gongzhuling, Jilin, China). 
Hemiechinus dauuricus (Sundevall), Satunin, 1907: 185. 
Erinaceus (Mesechinus) dauuricus (Sundevall), Ognev, 1951: 
8-14. 

Mesechinus dauuricus (Sundevall), Frost et al., 1991: 30. 


Hedgehog of genus Mesechinus (GLS-55.18 mm; Table 1). 


Length of ear similar to surrounding spines. Spines 21-23 mm 
in length, white for two-thirds of length, followed by black ring, 
light narrow ring, and black tip (Figure 3); premaxilla extending 
posteriorly to frontal (Figure 4); P3 triangular (equal-sided) in 
shape, similar to P? in size (Figure 4). 


Distribution: Widely distributed in eastern Inner Mongolia, 
Shaanxi, Ningxia, Heilongjiang, Jilin, and Liaoning, China; NE 


Mongolia; Transbaikalia and upper Amur Basin, Russia (Figure 1). 


Comments: Because M. miodon has been recognized 
previously as a synonym or subspecies of this species, the 
distribution boundary between these two species is unclear and 
the specimens from the southern-most distributions (especially 
in northwestern China Ningxia and Shaanxi) need to be 
carefully re-examined. 


Mesechinus hughi (Thomas, 1908) Hugh’s Forest Hedgehog 
Erinaceus hughi Thomas, 1908: 966. Type locality: Paochi 
(=Baoji), Shaanxi, China. 

Hemiechinus sylvaticus Ma, 1964: 31-36. Type locality: 
Qin-Shui District, Northern slope of Mt. Lishan, Shanxi, China. 
Hemiechinus dauuricus (Sundevall), Corbet, 1978: 15. 
Mesechinus hughi (Thomas), Frost et al., 1991: 30-31 
(including M. sylvaticus). 


Smallest species of Mesechinus (GLS-48.46 mm; Table 1). 


Ears not longer than surrounding spines; ventral pelage light 
brown; spines 19—21 mm in length, color pattern same as 
M. dauuricus (Figure 3); frontal relatively higher than parietal; 
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short spine on posterior palatal shelf moderately developed; 
epipterygoid process moderately developed; basisphenoid 
moderately inflated; nasal narrow, premaxilla extending 
posteriorly, frontal extending anteriorly, not meeting premaxilla, 
nasal and maxilla sharing short suture; P? triangular (nearly 
equal-sided) in shape, smaller than P2. 


Distribution: Southern Shaanxi, southern Shanxi, and 
northern Sichuan in China (Figure 1). 


Comments: We recognized H. sylvaticus (Ma, 1964) as a 
synonym of M. hughi. To date, its taxonomic status has 
not been appropriately evaluated as it is only known from its 
holotype. Ma (1964) described sylvaticus as a new species 
but did not examine the specimens of M. hughi. Here, the 
characters used to define sylvaticus, such as spine color 
pattern and presence of sagittal ridge, were observed on all 
specimens of M. hughi examined. Therefore, we recognized 
sylvaticus as a synonym of M. hughi. This species inhabits 
mountainous broad-leaved forest, distinct from M. dauuricus 
and M. miodon, and its overall dark color may be an adaptation 
to such environments. 


Mesechinus miodon (Thomas, 1908) Small-toothed Forest 
Hedgehog 

Erinaceus miodon Thomas, 1908: 965. Type locality: Yulinfu 
(2Yulin), Shaanxi, China. 

Erinaceus europaeus miodon Thomas, Allen, 1938: 47—54. 
Hemiechinus dauuricus (Sundevall), Corbet, 1978: 15. 
Hemiechinus dauuricus (Sundevall), Min & Lin, 1989: 4. 
Mesechinus dauuricus (Sundevall), Frost et al., 1991: 30. 
Mesechinus miodon (Thomas), Wang, 2003: 4. 


Large species of Mesechinus (GLS=54.10 mm; Table 
1) Ventral pelage pale white; spines 22-29 mm long, 
first two-thirds light brown (ivory white in other species), 
then broadly ringed blackish-brown, terminal 3-4 mm of 
spine light brown; parietal higher than frontal; well-developed 
spine on posterior palatal shelf; epipterygoid processes well 
developed; nasal narrow, premaxilla extending posteriorly, 
frontal extending anteriorly on dorsal side, touching premaxilla, 
nasal and maxilla without common suture; P? triangular (nearly 
equal-sided) in shape, smaller than P2. 


Distribution: Northern Shaanxi and eastern Ningxia, China 
(Figure 1). 


Comments: Named as Erinaceus miodon based on small P3 
(Thomas, 1908). Corbet (1978) considered it as a synonym 
of Hemiechinus dauuricus, which was subsequently followed 
by many researchers (e.g., Corbet, 1978, 1988; Corbet & Hill, 
1992). Hutterer (2005) assigned it as a synonym of M. hughi 
(perhaps following a comment in Frost et al. (1991)). However, 
both the skull size and shape were very distinct from M. hughi 
in our morphometric analyses (Figure 6). Furthermore, it could 
be distinguished from M. dauuricus based on the different 
color patterns on the spine (Figure 3) and smaller P3. The 
epipterygoid processes were also longer than those in M. 


dauuricus, though the sample size was small (Figure 4). We 
therefore recognized it as a valid species. 

The karyotype of animals from the type locality varied from 
2n=44—48 due of the presence of B chromosomes (Lin & Min, 
1989; Kong et al., 2016a). Kong et al. (2016b) reported 
both M. miodon and M. dauuricus as distributed in Yuling 
in northern-most Shaanxi; however, the distribution boundary 
remains unknown. 


Mesechinus wangi sp. nov. HE, JIANG, and AI 

Common names: Gaoligong Forest Hedgehog (im £& Di Jk 
7H, Gaoligong Linwei) or Wang's Forest Hedgehog (E KÂVE, 
Wangshi Linwei) 














Holotype: KIZ 022028 (field number: 201012001), adult 
female collected from Gaoligongshan National Nature Reserve 
(N24°50’, E98°45’), Baoshan, Yunnan, China, on 1 September 
2010 at an altitude of 2 215 m as... Alcohol-preserved and 
cleaned skull are deposited in KIZ, CAS. 


Paratypes: KIZ 027001 (field number: 0907003), KIZ 027002 
(field number: 0907001), KIZ 022027 (field number: 1102007), 
KIZ 034255 (field number: 201507001), and KIZ 034115 
(field number: GLGS 20160601) collected from Gaoligongshan 
National Nature Reserve, southwestern Yunnan, China from 
2003 to 2016 at elevations of 2 100 to 2 680 m. Except for KIZ 
034115, which is preserved in fluid, all other specimens are 
preserved as dried skins with cleaned skulls. The skull of KIZ 
027002 is broken and the skull of KIZ 034255 is missing. 


Etymology: Named in memory of Prof. Ying-Xiang Wang 
(1938-2016), head of the mammal research group at the 


Kunming Institute of Zoology, Chinese Academy of Sciences. 


He undertook extensive research on the taxonomy, phylogeny, 
zoogeography, and conservation of mammals and made 
distinguished contributions to Mammalogy in China (Jiang, 
2016). 


Diagnosis: Body size larger than M. hughi, but similar to M. 
dauuricus and M. miodon. Most spines (28096) lack white 
ring, in contrast to other Mesechinus species. Frontal higher 
than parietal, differing from M. dauuricus and M. miodon 
but similar to M. hughi. Spine on posterior palatal shelf 
short, similar to M. hughi but different from M. dauuricus 
and M. miodon.  Epipterygoid processes longer than that 
on M. dauuricus and M. hughi, but shorter than that on 
M. miodon. | Basisphenoid uninflated, distinct from other 
Mesechinus species with basisphenoids moderately inflated. 
Nasal (74.30 mm) broader than all other Mesechinus species 
(«3.00 mm). P? smaller than that of M. dauuricus, but similar 
to M. hughi and M. miodon. Supernumerary M4 consistently 
present after each M?, unique among all living hedgehogs 
(Figure 5). 


Description: Large Mesechinus species (HB=202.40 mm; 
GLS=54.75 mm; CBL=54.55 mm; Table 1). Absence of 
spineless area on scalp; length of ear equal to surrounding 
spines; ventral pelage dark brown; spines 22-25 mm long, 
most spines (28096) white for two-thirds of length and black 


for other one-third, small number of spines (<20%) white for 
two-thirds of length, then ringed in black, followed by narrow 
white ring, tip black; frontal higher than parietal; short spine 
("1 mm) present on posterior palatal shelf, extending only 
slightly posteriorly; epipterygoid processes extend labially (2-3 
mm); basisphenoids uninflated, two basisphenoids on both 
sides touch medially, excavated into shallow spherical fossa 
(namely, nasopharyngeal fossa; see Frost et al., 1991), breadth 
of nasopharyngeal fossa ~1.5 mm, breadth of nasal ~4.3 mm; 
premaxilla extending slightly posteriorly and frontal extending 
slightly anteriorly on dorsal side of skull, not meeting each 
other, nasal and maxilla sharing suture (~5—8 mm); jugal large, 
reaching lacrimal, lacrimal/maxilla suture unfused in adults; 
I? small, I? with two roots, larger than l?, P? rectangular 
shaped, similar to I? in size but larger than C', P? smaller 
than P2, triangular (nearly equal-sided) in shape, M? heavily 
reduced, hypocone and metacone absent, M^ single rooted, 
much smaller than M? (Figure 5). 


Measurement: Measurements for Mesechinus wangi sp. nov. 
(KIZ 027028, 027001, 027002, 022027 and 034225) are 
presented in Table 4. 


Comparisons: Mesechinus wangi sp. nov. can be 
characterized by many unique features within Mesechinus, 
including unique color pattern on spine,  uninflated 
basisphenoid, broad nasal, long common suture shared by 
nasal with maxilla, and presence of M^. It is similar to M. 
dauuricus and M. miodon in overall size (HB=202.40+26.10 
mm) but is obviously larger than M. hughi (HB=189.71+24.20 
mm; Table 1). Mesechinus wangi sp. nov. differs from 
M. dauuricus and M. miodon in relatively higher frontal than 
parietal in skull. 


Distribution: To date, this species is known only from three 
counties (Tengchong, Longling, and Longyang) of Baoshan in 
Yunnan, China, at elevations ranging from 2 200 m-2 680 m. 
The habitat is subtropical evergreen broad-leaved forest formed 
by a variety of vegetation, including Fagaceae, Lauraceae, 
Ericaceae, and Theaceae (Figure 8). 


Comments: Population size is currently unknown. However, 
the known distribution is extremely small and located only 
within the Gaoligong National Nature Reserve. The species 
hibernates from middle of October to the following early April. 


Key of four species of Mesechinus 
1. Ventral pelage dark brown; white for two-thirds and 
black for other one-third on most spines; greatest width 
of nasal>4.00 mm; nasal shares long common suture 
with maxilla; basisphenoid uninflated; supernumerary M^ 
present after M?, occurs only in Mt. Gaoligong, Yunnan, 
China..…...................…..Mesechinus wangi sp. nov. 


Ventral pelage light brown or white; nasal narrower 
than 3.00 mm; premaxilla meets frontal on dorsal 
side of skull (or nearly), nasal does not share 
suture with maxilla; basisphenoid moderately 
inflated; M^ not present; occurs outside of 
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2. Overall small; GLS<53.00 mm, LUTR<27.00 mm; frontal 
relatively higher than parietal; occurs in southern 
Shaanxi, southern Shanxi, and northern Sichuan, 


Overall large; GLS>53.00 mm, LUTR>27.00 mm; parietal 
relatively higher than frontal, spine on posterior palatal 
shelf well developed........................................3 


3. Spines 21-23 mm in length; tip of spines 
black, followed by narrow white ring; epipterygoid 
processes short; P? similar to P? in size; 
distributed in northern China, Mongolia, and 
RUSSIE D M. dauuricus 


Spines 22-29 mm in length; tip of spine light brown; 
epipterygoid processes well developed; P? obviously 
smaller than P2..................................M. miodon 


Table 4 External and craniodental measurements for type 
specimens of Mesechinus wangi sp. nov. 








: Holotype Paratypes 

Variable 

022028 022027 027001 027002 034255 
W 449.00 336.00 390.00 430.00 451.00 
HB 240.00 215.00 200.00 180.00 177.00 
TL 18.20 18.00 14.00 18.10 18.00 
HF 48.00 48.00 45.30 46.70 48.00 
EL 31.10 28.90 31.80 28.20 28.00 
GLS 55.10 55.60 53.70 54.60 — 
CBL 54.80 55.20 53.60 54.60 — 
CH 17.50 17.60 16.10 17.30 — 
BL 50.50 51.30 50.50 47.70 — 
PL 30.80 30.60 30.10 29.50 — 
ZMB 34.10 33.70 34.10 — — 
IOB 15.10 14.80 14.20 14.60 — 
MTW 26.20 26.20 25.30 24.70 — 
GWN 4.30 4.30 4.30 - — 
BM! 21.50 21.10 21.70 — - 
LUTR 29.10 28.70 27.10 26.70 — 
LBTR 25.30 24.20 25.20 24.70 — 


Abbreviations are explained in the Materials and Methods section. 
—: Not available. 
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ABSTRACT 


Karyotypes of four Chinese species of field mice 
of the genus Apodemus were examined, including 
Apodemus chevrieri (diploid chromosome number, 
2n-48, fundamental number of autosomal arms, 
FNa-56), A. draco (2n-48, FNa=48), A. ilex 


(2n=48, FNa=48), and A. latronum (2n=48, FNa=48). 


Karyotypes of A. chevrieri, A. draco, and A. ilex 
are reported here for the first time, providing useful 
information for their species taxonomy. Determining 
the karyotypes of all species of Apodemus in Asia, 
both in this and previous studies, provides a solid 
overview of the chromosome evolution and species 
differentiation of the genus in East Asia. In addition 
to allopatric speciation, chromosome rearrangements 
likely played an important role in the formation of the 
four Apodemus species groups as well as speciation 
within each group in East Asia. For example, increased 
centromeric heterochromatin in A. latronum may have 
contributed to the post-mating reproductive isolation 
from the A. draco-A. ilex-A. semotus clade. 


Keywords: Karyotype; Chromosome evolution; 
Speciation; Taxonomy; Field mice 


INTRODUCTION 


Field mice of the genus Apodemus are common murid species 
widely distributed in the Palearctic region through to the 
northern part of the Oriental region. The genus currently 
includes 20 species (Musser et al., 1996; Musser & Carlenton, 


348 Science Press 


2005), which have been characterized into three species 
groups based on morphological characters from detailed 
literature review (Musser et al., 1996): that is, Apodemus 
Group (A. agrarius, A. chevrieri, A. speciosus, A. peninsulae, 
A. latronum, A. draco, A. semotus, A. gurkha), Sylvaemus 
Group (A. sylvaticus, A. flavicollis, A. uralensis, A. mystacinus, 
A. fulvipectus, A. heremonensis, A. alpicola, A. arianus, A. 
hyrcanicus, A. ponticus, A. rusiges, A. wardi), and Argenteus 
Group (A. argenteus). The Apodemus Group and Argenteus 
Group consist of species distributed in East Asia, whereas 
species within the Sylvaemus Group are found in western 
Palearctic region. The A. agrarius species from the Apodemus 
Group is widely distributed in the Palearctic region from East 
Asia to Europe. Currently, however, there is still considerable 
taxonomic confusion regarding the species boundaries and 
identification of East Asian Apodemus species (Musser et al., 
1996), especially those distributed in China. 

Several phylogenetic studies using genetic approaches were 
conducted to reveal the species relationship and validity of 
the above-mentioned species groups (Filippucci et al., 2002; 
Liu et al., 2004; Michaux et al., 2002; Serizawa et al., 2000; 
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Suzuki et al., 2003, 2008). Suzuki et al. (2008) conducted 
comprehensive phylogenetic analyses based on mitochondrial 
and nuclear genes from most species of Apodemus and 
confirmed the distinct lineages of the three species groups, 
except for A. gurkha, which showed an independent lineage 
from the other species within the Apodemus Group. 

Concerning the evolutionary history of the genus Apodemus 
in East Asia, Suzuki et al. (2008) determined that the three 
species groups formed around 6 million years ago (Ma), with 
the Apodemus Group splitting into four ancestral species (A. 
agrarius/A. chevrieri, A. draco (and A. ilex)/A. semotus/A. 
latronum, A. peninsulae, and A. speciosus) around 5 Ma, and 
then splitting into the currently recognized species around 2 Ma. 
For these speciation events, Suzuki et al. (2008) assumed that 
allopatric speciation likely played an important role, followed by 
range expansion and distribution overlap. The original place 
for speciation event, however, has not been mentioned and 
unspecified. 

Chromosomal divergence is thought to play a role in 
reproductive isolation (e.g., King, 1993). Examination 
of karyotypes of species and populations is important to 
reconstruct allopatric and sympatric speciation events and 
clarify the historical changes in species distribution. Species 
differentiation among congeneric species also participates in 
cytological reproductive isolation (e.g., King, 1993). While the 
karyotypes of Apodemus species have been relatively well 
studied (e.g., Matsubara et al., 2004), information on species 
and populations in China is still limited. Clarification of species 
karyotypes is important for understanding the diversification 
of a genus. In this study, we examined the karyotypes of 
A. chevrieri, A. draco, A. ilex, and A. latronum based on 
specimens collected in China to help fill the gap in current 
knowledge. Even though the newly reported karyotypes were 
limited to conventional karyotypes, we expect they will be useful 
for the evaluation of species taxonomy and will provide an 
overview of chromosomal evolution and species differentiation. 
We also examined evolutionary history in consideration of the 
molecular and chromosomal divergences of Apodemus in East 
Asia. 


MATERIALS AND METHODS 


A total of 71 specimens from four Apodemus species (A. 
chevrieri, A. draco, A. ilex, and A. latronum) in China 
were examined. Species identification was made by careful 
examination of cranial characters following Musser et al. 
(1996), in addition to external characters and measurements. 
Apodemus ilex (mostly distributed in Yunnan, China) is often 
considered a synonym of A. draco (e.g., Musser & Carlenton, 
2005); however, molecular phylogeographic data suggest 
two species (e.g., Liu et al, 2012). In this study, we 
considered A. ilex as a separate species from A. draco, 
even though future study is expected to evaluate their 


taxonomic status and geographic distribution more accurately. 


Voucher specimens were deposited in the Key Laboratory of 
Conservation and Application in Biodiversity of South China, 
Guangzhou University, Guangzhou (GU), and the Marine 


College of Shandong University at Weihai (SUS). 


Examined specimens and collection localities are as follows: 
Apodemus chevrieri (n=11): Mt. Emei, Sichuan, GU 
MM3566 (male), 3593, 3594, 4478, 4480, 4484 (females), 
Wolong, Sichuan, SUS $1124, 81264, 81265 (males), S1107, 
$1236 (females); Apodemus draco (n-41) Mt. Emei, 
Sichuan, GU MM3545, 3563, 3564, 3568, 3569, 3570, 3585, 
3586, 3596, 3599, 4479, 4483, 4485 (males), 3551, 3565, 
3578, 3579, 3587, 3595, 4482 (females); Labahe, Tianquan, 
Ya'an, Sichuan, GU10073, 10076, 10077, 10094, 10107, 
10128 (males), 10074, 10108, 10110 (females); Kangding, 
Sichuan, GU10137, 10139, 10148 (males), 10135, 10147 
(females); Wolong, Sichuan, SUS S1140, $1257, S1266 
(males), 81108, S1180, S1245, S1246 (females); Apodemus 
ilex (n=9): Ailaoshan, Xinping, Yunnan, SUS S570, S649, 
$661, S663, S667, S683 (males), S651, S662, S684 (females); 
Apodemus latronum (n=10): Kangding, Sichuan, GU10134, 
10157 (males), 10136, 10140, 10145, 10151, 10153 (females), 
Wolong, Sichuan, SUS 81136, 81156 (males), S1134 (female). 

Cytological preparations were made from tail and/or lung 
tissue culture cells using the standard air-drying method 
described by Harada & Yosida (1978). C-band staining was 
accomplished as per Sumner (1972) for selected species and 
specimens. Terminology for chromosomes followed Levan et 
al. (1964): i.e., metacentric, submetacentric, subtelocentric, 
and acrocentric. Diploid chromosome number (2n) and 
fundamental number of autosomal arms (FNa) were calculated. 


RESULTS 


The karyotype of Apodemus chevrieri (Figure 1A) consisted 
of four small meta- or submetacentric pairs (nos. 1—4) and 
19 large-to-small acrocentric pairs (nos. 5—23) in autosomes, 
large acrocentric X chromosome, and small acrocentric Y 
chromosome. The 2n and FNa values were 48 and 54, 
respectively. 

The karyotype of Apodemus draco (Figure 1B) consisted 
of one small metacentric pair (no. 1) and large-to-small 
acrocentric pairs (nos. 2—23) in autosomes, large acrocentric 
X chromosome, and small acrocentric Y chromosome. The 2n 
and FNa values were 48 and 48, respectively. 

The karyotype of Apodemus ilex (Figure 1C) consisted 
of one small metacentric pair (no. 1) and large-to-small 
acrocentric pairs (nos. 2—23) in autosomes, large acrocentric 
X chromosome, and small acrocentric Y chromosome. The 2n 
and FNa values were 48 and 48, respectively. 

The karyotype of Apodemus latronum (Figure 1D) consisted 
of one small submetacentric (no. 1) and 22 large-to-small 
(nos. 2-23) acrocentric pairs in autosomes, large acrocentric X 
chromosome, and small acrocentric Y chromosome. In several 
acrocentric pairs, the centromeric region was well developed 
due to the constitutive heterochromatins, which were well 
stained following C-band staining (Figure 1E, nos. 2-9). As we 
could not find clear short arms for those pairs, we considered 
those pairs to be acrocentric. The 2n and FNa values were 48 
and 48, respectively. 
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Figure 1 Karyotypes of Apodemus species from China 
Conventional karyotypes of A. chevrieri (A, GU MM3593), A. draco (B, 
M10077), A. ilex (C, SUS S649), and A. latronum (D, GU10134), as well 


as the C-band karyotype of A. latronum (E, GU10134). 


DISCUSSION 


We analyzed the karyotypes of four Apodemus species 
from China. Previous karyotypic data from this genus are 
summarized in Table 1, together with our results from this study. 

The karyotype of A. chevrieri is reported here for the 
first time, and was characterized by four small metacentric 
pairs (2n248, FNa-54). Apodemus chevrieri is restricted to 
southwestern China and based on mitochondrial and nuclear 
gene phylogenetic studies is thought to be a sister or in-group 
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species of the widely distributed A. agrarius (Liu et al., 
2004; Suzuki et al., 2003, 2008). Although the karyotype 
of A. agrarius is polymorphic and possesses 3-5 biarmed 
metacentric autosome pairs (2248, FNa=52—56, excluding 
the B chromosome; Boeskorov et al., 1995; Britton-Davidian 
et al., 1991; Chassovnikarova et al., 2009; Chernukha et al., 
1986; Kang & Koh, 1976; Kartavtseva, 1994; Kartavtseva & 
Pavlenko, 2000; Kefelioglu et al., 2003; Koh, 1987, 1988, 1989; 
Král, 1970, 1972; Matsubara et al., 2004; Shbulatova et al., 
1991; Soldatović et al., 1969, 1975; Tsuchiya, 1979; Vujošević 
et al., 1984; Wang et al., 1993; Yigit et al., 2000), the karyotype 
with four metacentric pairs (2n248, FNa-54) is regarded 
as the standard karyotype for A. agrarius (see Kartavtseva 
& Pavlenko, 2000). Therefore, we suggest that there are 
no clear differences in the conventional karyotypes between 
A. chevrieri and A. agrarius; however, further study using 
differential staining of chromosome arms is expected to clarify 
any minor differences and rearrangement of chromosome arms 
between A. chevrieri and polymorphic A. agrarius, and thus 
help reevaluate their taxonomic status. 

The karyotypes of A. draco and A. ilex are reported in 
this study for the first time as correct species identification, 
with both characterized by one small metacentric pair (20-48, 
FNa-48), similar to that of A. semotus in Taiwan, China 
(Matsubara et al., 2004; Tsuchiya, 1979). While Chen 
et al. (1996) reported karyotypes of A. draco as 2n-48, 
FNa=46 and A. peninsulae as 2n=48, FNa=48 from Yunnan 
Province, China, these two karyotypes were possibly reported 
based on erroneous identification. We suggest that the 
former specimens collected from Kunming were A. peninsulae, 
whereas the latter specimens collected from Jianchuan were 
A. ilex. This interpretation of misidentification by Chen et 
al. (1996) would be congruent with the distribution of A. 
draco (currently A. ilex) in Kunming and Jianchuan and A. 
peninsulae in Kunming but not Jianchuan (Zhang, 1997); and 
that these two species have been considered superficially 
similar in morphologies and often misidentified before the 
careful taxonomic revision by Musser et al. (1996). 

The karyotype of specimens of "A. draco" by Chen et al. 
(1996), and herewith interpreting to represent A. peninsulae 
showed no differences with the reported A. peninsulae 
karyotype and had only acrocentric chromosomes (2n=48, 
FNa=46; Hayata, 1973; Kartavtseva et al., 2000; Koh, 1986, 
1988; Wang et al, 2000). The karyotype of the latter 
specimens correctly representing A. ilex was very similar 
to the karyotype for A. ilex from Yunnan, as well as A. 
draco from Sichuan in this study (2n248, FNa-48) and A. 
semotus from Taiwan, China (2n=48, FNa-48; Matsubara et al., 
2004; Tsuchiya, 1979) characterized by one small metacentric 
pair. Although the current study was limited to conventional 
karyotypes, we report here on the karyotypes of A. draco and 
A. ilex for the first time and provide updated information on 
the karyotype of A. peninsulae. These data are important 
for further study on species taxonomy and identification of the 
genus Apodemus in East Asia. 


Table 1 Karyotypes of field mice of the genus Apodemus examined in this study and reported in previous studies 





Species Locality 2n FNa M/SM ST A X Y B Reference 
A. chevrieri Sichuan, China 48 54 4 0 20 A A - This study 
A. agrarius Shandong, China 48 54 4 0 19 A A - Wang et al. (1993) 
Taiwan, China 48 56 5 0 18 A A - Tsuchiya (1979) 
Korea 48 54 4 0 19 A A - Kang & Koh (1976), Koh (1987, 1988, 1989), 
Matsubara et al. (2004) 
Primorye 48 52 3 0 20 A A - Chernukha et al. (1986) 
Primorye 48 52-54 3-4 0 19-20 A A 0-1 Kartavtseva & Pavlenko (2000) 
Amur 48 52 3 0 20 A A - Kartavtseva & Pavlenko (2000) 
Khasan 48 54 4 0 19 A A - Boeskorov et al. (1995) 
Khabarovsk 48 52-54 3-4 0 19-20 A A 0-1 Chernukha et al. (1986), Kartavtseva (1994), 
Kartavtseva & Pavlenko (2000) 
Siberia 48 52-54 3-4 0 19-20 A A - Boeskorov et al. (1995), Kartavtseva & Pavlenko (2000) 
Altai 48 52 3 0 20 A A - Chernukha et al. (1986) 
Altai 48 54 4 0 19 A A - Kartavtseva & Pavlenko (2000) 
Moskow oblast 48 52 3 0 20 A A - Chernukha et al. (1986) 
Chechen-Ingush 48 52 3 0 20 A A - Chernukha et al. (1986) 
Krasnodar 48 52 3 0 20 A A - Chernukha et al. (1986) 
Ukraine 48 54 4 0 19 A A - Kartavtseva & Pavlenko (2000) 
Moldova 48 52-54 3-4 0 19-20 A A - Kartavtseva & Pavlenko (2000) 
Azerbaijan 48 54 4 0 19 A A - Shbulatova et al. (1991) 
Czechoslovakia 48 54 4 0 19 A A - Kral (1970) (1972) 
Poland 48 54 4 0 19 A A - Kral (1970) 
Yugoslavia 48 54 4 0 19 A A - VujoSevié et al. (1984) 
Yugoslavia 48 52-54 3-4 0 19-20 A A - Soldatovié et al. (1969, 1975) 
Bulgaria 48 52-54 3-4 0 19-20 A A 0-1 Chassovnikarova et al. (2009) 
Greece 48 54 4 0 19 A A - Britton-Davidian et al. (1991) 
Turkey 48 54 4 0 19 A A - Kefelioglu et al. (2003) 
Turkey 48 56 5 0 18 A A - Yigit et al. (2000) 
A. draco Sichuan, China 48 48 1 0 22 A A - This study 
A. ilex Yunnan, China 48 48 1 0 22 A A - This study 
Yunnan, China 48 48 1 0 22 A A - Chen et al. (1996) as “A. peninsulae” 
A. latronum Sichuan, China 48 48 1 0 22 A A - This study 
Yunnan, China 48 66 8 2 13 A ? - Chen et al. (1996) 
A. semotus Taiwan, China 48 48 1 0 22 A ? - Matsubara et al. (2004) 
A. peninsulae Yunnan, China 48 46 0 0 23 A A - Chen et al. (1996) as “A. draco” 
NE China 48 46 0 0 23 A A 0-14 Wang etal. (2000) 
Korea 48 46 0 0 23 A A 6-1 Koh (1986, 1988) 
Russia 48 46 0 0 23 A A 0-6  Kartavtseva et al. (2000) 
Hokkaido, Japan 48 46 0 0 23 A A 0-13 Hayata (1973) 
A. speciosus Japan 46-48 54 4-3 1 17-19 A A - Tsuchiya (1974) 
Japan 46-48 54 5-4 0 17-19 A A - Saitoh & Obara (1986) 
A. argenteus Japan 46 50 2 0 20 SM A 0-1 Yoshida et al. (1975), Obara & Sasaki (1997) 
A. gurkha Nepal 48 50 2 0 21 A ? - Matsubara et al. (2004). 
Nepal 48 62-64 4-3 5 14-15 A A - Gemmeke & Niethammer (1982) 
Sylvaemus Group 
A. sylvaticus 48 46 0 0 23 A A - Zima & Král (1984), Orlov et al. (1996), 
KryStufek & Vohralík (2009) 
A. flavicollis 48 46 0 0 23 A A 1-3  Zima & Kral (1984), Orlov et al. (1996), 
KryStufek & Vohralík (2009) 
A. microps 48 46 0 0 23 A A - Zima & Král (1984), Reutter et al. (2001) 
A. alpicola 48 46 0 0 23 A A - Reutter et al. (2001) 
A. witherbyi 48 46 0 0 23 A A - Orlov et al. (1996), KryStufek & Vohralík (2009) 
A. uralensis 48 46 0 0 23 A A - Orlov et al. (1996), KryStufek & Vohralik (2009) 
A. ponticus 48 46 0 0 23 A A - Orlov et al. (1996) 
A. pallipes 48 46 0 0 23 A A - Gemmeke & Niethammer (1982) 
A. epimelas 48 48-50 1-2 0 21-22 A A 0-1 Belcheva et al. (1988), Zima & Kral (1984) 
A. mystacinus 48 50 2 0 21 A A - KryStufek & Vohralik (2009) 


Diploid and sex chromosomes were classified into metacentric (M), submetacentric (SM), subtelocentric (ST), and acrocentric (A), and a "?" 


indicate the Y chromosome was too small to be confirmed. 2n and FNa, excluding the B chromosome. —: Not available. 
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Figure 2 Hypothesized diversification process with allopatric distribution and chromosome changes in the genus Apodemus in 


East Asia 


Phylogenetic relationships among species followed the molecular phylogeny of Suzuki et al. (2008). Some chromosome rearrangements referred to Matsubara 


et al. (2004). Arrowheads and closed circles indicate possible chromosome rearrangements and allopatric speciation, which resulting reproductive isolation. ?: 


Indicate the hypothetical origination of the clade/lineage. 


The karyotype of A. latronum was 2n=48 and FNa=48, with 
one small biarmed pair. This chromosome complement was 
similar to that of A. draco, A. ilex, and A. semotus, but the 
karyotype differed by having centromeric heterochromatin in 
many acrocentric pairs. Similar centromeric heterochromatin 
has been found in previous study on the karyotype of A. 
latronum from Yunnan Province (Chen et al., 1996). Chen et 
al. (1996) stated that the centromeric heterochromatin formed 
short arms and thus considered the A. /atronum karyotype to 
be 2n=48, FNa=66. Although we did not analyze the G-band 
karyotype of A. latronum, based on the C-band karyotype we 
found no considerable differences between our A. latronum 
karyotype (2n=48, FNa=48) and that of Chen et al. (1996) 
(2n=48, FNa=66), despite different FNa values due to the 
interpretation of centromeric heterochromatin. 

We studied the karyotypes of all Apodemus species in East 
Asia and provided a solid overview of chromosome evolution 
and species differentiation of the genus within East Asia. The 
chromosome rearrangements in East Asian Apodemus were 
congruent with the species divergence pattern proposed in 


previous molecular study (Suzuki et al., 2008). Suzuki et al. 


(2008) recognized four groups as the major DNA phylogenetic 
clades of the East Asian Apodemus subgeneric group: (1) 
A. agrarius-A. chevrieri (=agrarius species group), (2) A. 
draco—A. ilex-A. semotus—A. latronum (=draco species group), 
(3) A. peninsulae, and (4) A. speciosus. Suzuki et al. (2008) 
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stated that these four groups radiated 6 Ma in response to 
global environmental changes among allopatric populations. 
Our present study clarified that these four DNA phylogenetic 
species groups were distinct, with different karyotypes: 21-48, 
FNa-54 for the agrarius group (A. agrarius, A. chevrieri); 
2n-48, FNa-48 for the draco group (A. draco, A. ilex, A. 
semotus, A. latronum); 2n=48, FNa=46 for A. peninsulae; and 
2n-46/48, FNa-54 for A. speciosus (Tsuchiya, 1974; Saitoh 
& Obara, 1986. We suggest that these major chromosome 
rearrangements among clades played an important role in 
clade formation through post-mating reproductive isolation, in 
addition to allopatric distribution. 

After the radiation into four groups, further speciation events 
are thought to have occurred within the draco and agrarius 
groups around 2 Ma (Suzuki et al., 2008). In the draco 
group, speciation likely occurred through allopatric speciation 
due to partitioning of the distribution range in developping 
geographic barriers, such as among A. ilex (Yunnan), A. draco 
(other areas in mainland China), and A. semotus (Taiwan, 
China), with minor chromosome rearrangements unlikely to 
have contributed to the speciation events of these three 
allopatric species (Figure 2). On the other hand, the current 
distribution range between A. /atronum and A. draco and 
between A. latronum and A. ilex overlap (e.g., Musser et al., 
1996). This suggests that A. /atronum, which is distributed 
in the western provinces of Sichuan, Yunnan, Xizang, and 


Qinghai, as well as northern Myanmar (Musser & Carlenton, 
2005), was not derived through allopatric speciation among 
the draco group. We propose that speciation of A. latronum 
from the A. draco-A. ilex-A. semotus clade may have occurred 
as sympatric speciation, where chromosome rearrangements 
contributed to form post-mating reproductive isolation at the 
cytological level. The increased centromeric heterochromatin 
found in A. /atronum also influenced post-mating reproductive 
isolation from the A. draco-A. ilex-A. semotus clade, which 
lacked heterochromatin increase (Figure 2). On the other hand, 
A. agrarius and A. chevrieri in the agrarius group exhibit slight 
overlap in their current distribution ranges (Musser et al., 1996); 
and these two species may have undergone speciation by 
allopatric distribution, with subsequent expansion and overlap 
of their distribution ranges, as discussed by Suzuki et al. (2008). 
The speciation of A. chevrieri from A. agrarius is, therefore, 
suggested to have been accompanied by allopatric speciation 
events, and this evolutionary story may explain the lack of major 
karyotypic differences between the two species. 

In addition, extensive geographical divergences within the 
species have been reported for morphological and genetic 
traits in East Asian Apodemus species: e.g., A. chevrieri 
(Yue et al., 2012), A. agrarius (Sakka et al., 2010), A. draco 
(Fan et al., 2012; Kaneko, 2010, 2012, 2015; Sakka et al., 
2010), A. ilex (Kaneko, 2010, 2012, 2015; Liu et al., 2012), A. 
latronum (Kaneko, 2010, 2012, 2015; Li & Liu, 2014; Sakka 
et al., 2010), A. semotus (Hsu et al., 2001), A. peninsulae 
(Kaneko, 2010, 2012, 2015; Sakka et al., 2010; Serizawa et 
al., 2002), A. speciosus (Kageyama et al., 2009; Shintaku et 
al., 2012; Shintaku & Motokawa, 2016; Suzuki et al., 2004; 
Tomozawa et al., 2014; Tomozawa & Suzuki, 2008), and A. 
argenteus (Suzuki et al., 2004). These complex patterns 
are thought to have formed through geographic isolation and 
genetic exchange (e.g., A. speciosus between Robertsonian 
chromosome races; Shintaku & Motokawa, 2016; Suzuki et al., 
2004; Tomozawa & Suzuki, 2008) after the formation of each 
species. More comprehensive analyses using morphology, 
chromosomes, and DNA markers are expected to clarify the 
complex evolutionary history of the Apodemus genus in East 
Asia. The present study elucidated the evolutionary pattern of 
the Apodemus genus in East Asia with reference to the major 
chromosome rearrangements at the among-species level. 
Future study of major and minor chromosome rearrangements 
at the within-species level using various chromosome arm 
staining techniques is expected. The genus Apodemus may 
be considered a good wild animal model to understand the 
roles of reproductive isolation by allopatric distribution and 
chromosome rearrangement during speciation events. 


COMPETING INTERESTS 
The authors declare that they have no competing interests. 
AUTHORS’ CONTRIBUTIONS 


M.M. and Y.C.L. designed the study. M.M., Y.W., M.H., Y.S., X.L.J, and Y.L 
collected specimens. M.M. made species identification. M.M., M.H., Y.W., 


and Y.C.L analyzed karyotypes. M.M. and Y.S. made literature surveys. M.M. 
wrote the manuscript. Y.W., M.H., Y.C.L. revised the manuscript. All authors 
read and approved the final manuscript. 


REFERENCES 


Belcheva RG, Topaschka-Ancheva MN, Atanassov N. 1988. Karyological 
studies of five species of mammals from Bulgarian fauna. Comptes rendus 
de l'Academie bulgare des Sciences, 42: 125—128. 

Boeskorov GG, Kartavtseva IV, Zagorodniuk IV, Belianin AN, Liapunova 
EA. 1995. Nucleolus organizer regions and B-chromosomes of field mice 
(Mammalia, Rodentia, Apodemus). Genetika, 31(2): 185—192. (in Russian) 
Britton-Davidian J, Vahdati M, Benmedhi F, Gros P, Nancé V, Croset 
H, Guerassimov S, Triantaphyllidis C. 1991. Genetic differentiation in 
four species of Apodemus from southern Europe: A. sylvaticus, A. 
flavicollis, A. agrarius and A. mystacinus (Muridae, Rodentia). Zeitschrift für 
Sáugetierkunde, 56: 25-33. 

Chassovnikarova T, Atanasov N, Dimitrov H. 2009. Chromosome 
polymorphism in Bulgarian populations of the striped field mouse (Apodemus 
agrarius Pallas 1771). Comparative Cytogenetics, 3(1): 1-9. 

Chen ZP, Liu RQ, Li CY, Wang YX. 1996. Studies on the chromosomes 
of three species of wood mice. Zoological Research, 17(3): 347—352. (in 
Chinese) 

Chernukha YG, Evdokimova OA, Chechovich AV. 1986. Results of karyologic 
and immunobiological studies of the striped field mouse (Apodemus 
agrarius) from different areas of its range. Zoological Journal, 65: 471—475. 
(in Russian) 

Fan ZX, Liu SY, Liu Y, Liao LH, Zhang XY, Yue BS. 2012. Phylogeography of 
the South China field mouse (Apodemus draco) on the southeastern Tibetan 
Plateau reveals high genetic diversity and glacial refugia. PLoS One, 7(5): 
638184. 

Filippucci MG, Macholán M, Michaux JR. 2002. Genetic variation and 
evolution in the genus Apodemus (Muridae: Rodentia). Biological Journal 
of the Linnean Society, 75(3): 395—419. 

Gemmeke H, Niethammer J. 1982. Zur Characterisierung der Waldmáuse 
(Apodemus) Nepals. Zeitschrift für Sáugetierkunde, 47: 33-38. 

Harada M, Yosida TH. 1978. Karyological study of four Japanese Myotis bats 
(Chiroptera, Mammalia). Chromosoma, 65(3): 283-291. 

Hayata |. 1973. Chromosomal polymorphism caused by supernumerary 
chromosomes in the field mouse, Apodemus giliacus. Chromosoma, 42(4): 
403-414. 

Hsu FH, Lin FJ, Lin YS. 2001. Phylogeographic structure of the Formosan 
wood mouse, Apodemus semotus Thomas. Zoological Studies, 40(2): 
91-102. 

Kageyama M, Motokawa M, Hikida T. 2009. Geographic variation in 
morphological traits of the large Japanese field mouse, Apodemus speciosus 
(Rodentia, Muridae), from the Izu Island group, Japan. Zoological Science, 
26(4): 266—276. 

Kaneko Y. 2010. Identification of Apodemus peninsula, draco and A. 
latronum in China, Korea, and Myanmar by cranial measurements. Mammal 
Study, 35(1): 31—55. 

Kaneko Y. 2012. Horizontal and elevational distributions of Apodemus 
peninsulae, A. draco and A. latronum. Mammal Study, 37(3): 183—204. 


Zoological Research 39(5): 348-355, 2018 353 


Kaneko Y. 2015. Latitudinal geographical variation of external and cranial 
measurements in Apodemus peninsulae, A. draco, and A. latronum. 
Mammal Study, 40(3): 143—165. 

Kang YS, Koh HS. 1976. Karyotype studies on three species of the family 
Muridae (Mammalia; Rodentia) in Korea. Korean Journal of Zoology, 19(3): 
101-112. 

Kartavtseva IV. 1994. A description of the B chromosomes in the karyotype 
of the field mouse Apodemus agrarius. Tsitologiia i Genetika, 28(2): 96. 
Kartavtseva IV, Pavlenko MV. 2000. Chromosome variation in the striped 
field mouse Apodemus agrarius (Rodentia: Muridae). Russian Journal of 
Genetics, 36(2): 162—174. 

Kartavtseva IV, Roslik GV, Pavlenko MV, Amachaeva EY, Sawaguchi S, 
Obara Y. 2000. The B-chromosome system of the Korean field mouse 
Apodemus peninsulae in the Russian Far East. Chromosome Science, 4: 
21—29. 

Kefelioglu H, Tez C, Gündüz |. 2003. The taxonomy and distribution of 
Apodemus agrarius (Pallas, 1771) (Mammalia: Rodentia) in the European 
part of Turkey. Turkish Journal of Zoology, 27: 141—146. 

King M. 1993. Species Evolution: The Role of Chromosome Change. 
Cambridge: Cambridge University Press. 

Koh HS. 1986. Systematic studies of Korean rodents: II. A chromosome 
analysis in Korean field mice, Apodemus peninsulae peninsulae Thomas 
(Muridae, Rodentia), from Mungyong, with the comparison of morphometric 
characters of these Korean field mice to symptric striped field mice, A. 
agrarius coreaee Thomas. Korean Journal of Systematic Zoology, 2: 1—10. 
Koh HS. 1987. Systematic studies of Korea rodents: III. Morphometric and 
chromosomal analyses of striped field mice, Apodemus agrarius chemuensis 
Jones and Johnson, from Jeju-Do. Korean Journal of Systematic Zoology, 3: 
24-40. 

Koh HS. 1988. Systematic studies of Korean rodents: IV. Morphometric and 
chromosomal analyses of two species of the genus Apodemus (Muridae). 
Korean Journal of Systematic Zoology, 4: 103—120. 

Koh HS. 1989. Systematic studies of Korean rodents. V. Morphometric 
and chromosomal analyses on island populations of striped field mice 
(Apodemus agrarius coreae) in southwestern coasts of the Korean Peninsula. 
Korean Journal of Systematic Zoology, 5: 1—12. 

Král B. 1970. Chromosome studies in two subgenera of the genus 
Apodemus. Zoologicke Listy, 19(2): 119—134. 

Král B. 1972. Chromosome characteristics of Muridae and Microtidae 
from Czechoslovakia. Acta Scientiarm Naturalium Academiae Scientiarum 
Bohemoslovacae Brno, 6(12): 1—78. 

KryStufek B, Vohralík V. 2009. Mammals of Turkey and Cyprus. Rodentia II: 
Cricetinae, Muridae, Spalacidae, Calomyscidae, Capromyidae, Hystricidae, 
Castoridae. Koper: Zgodovinsko drustvo za juzno Primorsko. 

Levan A, Fredga K, Sandberg AA. 1964. Nomenclature for centromeric 
position on chromosomes. Hereditas, 52(2): 201—220. 

Li S, Liu SY. 2014. Geographic variation of the large-eared field mouse 
(Apodemus latronum Thomas, 1911) (Rodentia: Muridae) with one new 
subspecies description verified via cranial morphometric variables and 
pelage characteristics. Zoological Studies, 53: 23. 

Liu Q, Chen P, He K, Kilpatrick CW, Liu SY, Yu FH, Jiang XL. 2012. 
Phylogeographic study of Apodemus ilex (Rodentia: Muridae) in Southwest 


354 www.zoores.ac.cn 


China. PLoS One, 7(2): e31453. 

Liu XM, Wei FW, Li M, Jiang XL, Feng ZJ, Hu JC. 2004. Molecular phylogeny 
and taxonomy of wood mice (genus Apodemus Kaup, 1829) based on 
complete mtDNA cytochrome b sequences, with emphasis on Chinese 
species. Molecular Phylogenetics and Evolution, 33(1): 1—15. 

Matsubara K, Nishida-Umehara C, Tsuchiya K, Nukaya D, Matsuda Y. 
2004. Karyotypic evolution of Apodemus (Muridae, Rodentia) inferred from 
comparative FISH analyses. Chromosome Research, 12(4): 383—395. 
Michaux JR, Chevret P, Filippucci MG, Macholán M. 2002. Phylogeny of the 
genus Apodemus with a special emphasis on the subgenus Sylvaemus using 
the nuclear IRBP gene and two mitochondrial markers: cytochrome b and 
128 rRNA. Molecular Phylogenetics and Evolution, 23(2): 123—136. 

Musser GG, Brothers EM, Carleton MD, Hutterer R. 1996. Taxonomy and 
distributional records of Oriental and European Apodemus, with a review of 
the Apodemus-Sylvaemus problem. Bonner Zoologische Beitráge, 46(1—4): 
143-190. 

Musser GG, Carlenton MD. 2005. Superfamily muroidea. /n: Wilson DE, 
Reeder DM. Mammal Species of the World: A Taxonomic and Geographic 
Reference. 3rd ed. Baltimore: The Johns Hopkins University Press, 
894—1531. 

Obara Y, Sasaki S. 1997. Fluorescent approaches on the origin of B 
chromosomes of Apodemus argenteus hokkaidi. Chromosome Science, 
1(1): 1-5. 

Orlov VN, Bulatova NS, Nadjafova RS, Kozlovsky Al. 1996. Evolutionary 
classification of European wood mice of the subgenus Sylvaemus based 
on allozyme and chromosome data. Bonner Zoologische Beitráge, 46: 
191—202. 

Reutter BA, Nová P, Vogel P, Zima J. 2001. Karyotypic variation between 
wood mouse species: banded chromosomes of Apodemus alpicola and A. 
microps. Acta Theriologica, 46(4): 353—362. 

Saitoh M, Obara Y. 1986. Chromosome banding patterns in five intraspecific 
taxa of the large Japanese field mouse, Apodemus speciosus. Zoological 
Science, 3: 785—792. 

Sakka H, Quéré JP, Kartavtseva |, Pavlenko M, Chelomina G, Atopkin 
D, Bogdanov A, Michaux J. 2010. Comparative phylogeography of four 
Apodemus species (Mammalia: Rodentia) in the Asian Far East: evidence 
of Quaternary climatic changes in their genetic structure. Biological Journal 
of the Linnean Society, 100(4): 797—821. 

Serizawa K, Suzuki H, Tsuchiya K. 2000. A phylogenetic view on species 
radiation in Apodemus inferred from variation of nuclear and mitochondrial 
genes. Biochemical Genetics, 38(1—2): 27—40. 

Serizawa K, Suzuki H, lwasa MA, Tsuchiya K, Pavlenko MV, Kartavtseva 
IV, Chelomina GN, Dokuchaev NE, Han SH. 2002. A spatial aspect on 
mitochondrial DNA genealogy in Apodemus peninsulae from East Asia. 
Biochemical Genetics, 40(5—6): 149—161. 

Shbulatova N, Nadjafova RS, Kozlovsky I. 1991. Cytotaxonomic analysis of 
species of the genera Mus, Apodemus and Rattus in Azerbaijan. Journal of 
Zoological Systematics and Evolutionary Research, 29(2): 139—153. 
Shintaku Y, Kageyama M, Motokawa M. 2012. Morphological variation in 
external traits of the large Japanese field mouse, Apodemus speciosus. 
Mammal Study, 37(2): 113—126. 

Shintaku Y, Motokawa M. 2016. Geographic variation in skull morphology of 


the large Japanese field mice, Apodemus speciosus (Rodentia: Muridae) 
revealed by geometric morphometric analysis. Zoological Science, 33(2): 
132-145. 

Soldatović B, Dulić B, Savić i Desanka Rimsa I. 1969. Chromosomen zweier 
arten der Gattung Apodemus (A. agrarius und A. mystacinus-Mammaiia, 
Rodentia) aus Jugoslawien. Arhiv bioloskih nauka Beograd, 21: 27-32. 


Soldatović B, Savić |, Seth P, Reichenstein H, Tolksdorf M. 1975. 


Comparative karyological study of the genus Apodemus (Kaup 1829). Acta 
Veterinaria, 25(1): 1—10. 

Sumner AT. 1972. A simple technique for demonstrating centromeric 
heterochromatin. Experimental Cell Research, 75(1): 304—306. 


Suzuki H, Sato JJ, Tsuchiya K, Luo J, Zhang YP, Wang YX, Jiang XL. 
2003. Molecular phylogeny of wood mice (Apodemus, Muridae) in East Asia. 


Biological Journal of the Linnean Society, 80(3): 469—481. 


Suzuki H, Yasuda SP, Sakaizumi M, Wakana S, Motokawa M, Tsuchiya K. 


2004. Differential geographic patterns of mitochondrial DNA variation in two 
sympatric species of Japanese wood mice, Apodemus speciosus and A. 
argenteus. Genes & Genetic Systems, 79(3): 165—176. 


Suzuki H, Filippucci MG, Chelomina GN, Sato JJ, Serizawa K, Nevo E. 


2008. A biogeographic view of Apodemus in Asia and Europe inferred from 
nuclear and mitochondrial gene sequences. Biochemical Genetics, 46(5—6): 
329-346. 

Tomozawa M, Suzuki H. 2008. A trend of central versus peripheral 
structuring in mitochondrial and nuclear gene sequences of the Japanese 
wood mouse, Apodemus speciosus. Zoological Science, 25(3): 273—285. 
Tomozawa M, Nunome M, Suzuki H, Ono H. 2014. Effect of founding 
events on coat colour polymorphism of Apodemus speciosus (Rodentia: 
Muridae) on the Izu Islands. Biological Journal of the Linnean Society, 
113(2): 522-535. 


Tsuchiya K. 1974. Cytological and biochemical studies of Apodemus 
speciosus group in Japan. Journal of the Mammalogical Society of Japan, 
6(2): 67-87. (in Japanese) 

Tsuchiya K. 1979. Notes on breeding of wood mouse groups for laboratory 
animal. Report of Hokkaido Institute of Public Health, 29: 102—106. (in 
Japanese) 

Vujošević M, Rimsa D, Zivcovié S. 1984. Patterns of G- and C-bands 
distribution on chromosomes of three Apodemus species. Zeitschrift für 
Sáugetierkunde, 49: 234—238. 

Wang JX, Zhao XF, Wang XM. 1993. Studies of chromosome of striped field 
mouse Apodemus agrarius pallidior (Rodentia). Acta Theriologica Sinica, 
13(4): 283-287. (in Chinese) 

Wang JX, Zhao XF, Qi HY, Koh HS, Zhang L, Guan ZX, Wang CH. 
2000. Karyotypes and B chromosomes of Apodemus peninsulae (Rodentia, 
Mammalia). Acta Theriologica Sinica, 20(4): 289—296. 

Yiğit N, Verimli R, Sözen M, Çolak E, Özkurt S. 2000. The karyotype of 
Apodemus agrarius (Pallas, 1771) (Mammalia: Rodentia) in Turkey. Zoology 
in the Middle East, 20(1): 21—23. 

Yoshida MC, Sasaki M, Oshimura M. 1975. Karyotype and heterochromatin 
pattern of the field mouse, Apodemus argenteus Temminck. Genetica, 45(4): 
397-403. 

Yue H, Fan ZX, Liu SY, Liu Y, Song ZB, Zhang XY. 2012. A mitogenome of the 
Chevrier's field mouse (Apodemus chevrieri) and genetic variations inferred 
from the cytochrome b gene. DNA and Cell Biology, 31(4): 460—469. 

Zhang RZ. 1997. Distribution of Mammalian Species in China. Beijing: China 
Forestry Publishing House. (in Chinese) 

Zima J, Kral B. 1984. Karyotypes of European mammals Il. Acta Scientarium 
Naturalium Academiae Scientiarum Bohemoslovacae Brno, 18(8): 1—62. 


Zoological Research 39(5): 348-355, 2018 355 


ZOOLOGICAL RESEARCH 


Species identification of crested gibbons (Nomascus) in 
captivity in China using karyotyping- and PCR-based 


approaches 


Wen-Hui Nie’, Jin-Huan Wang', Wei-Ting Su’, Yu Hu!, Shui-Wang He!, Xue-Long Jiang’, Kai Het?” 


! Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming Yunnan 650223, China 


2 Kyoto University Museum, Kyoto University, Kyoto 606-8417, Japan 


ABSTRACT 


Gibbons and siamangs (Hylobatidae) are well-known for 
their rapid chromosomal evolution, which has resulted 
in high speciation rate within the family. On the other 
hand, distinct karyotypes do not prevent speciation, 
allowing interbreeding between individuals in captivity, 
and the unwanted hybrids are ethically problematic as all 


gibbon species are endangered or critically endangered. 


Thus, accurate species identification is crucial for captive 
breeding, particularly in China where studbooks are 
unavailable. Identification based on external morphology 
is difficult, especially for hybrids, because species 
are usually similar in appearance. In this study, we 
employed G-banding karyotyping and fluorescence in 
situ hybridization (FISH) as well as a PCR-based 
approach to examine karyotypic characteristics and 
identify crested gibbons of the genus Nomascus from 
zoos and nature reserves in China. We characterized 
and identified five karyotypes from 21 individuals of 
Nomascus. Using karyotypes and mitochondrial and 
nuclear genes, we identified three purebred species 
and three hybrids, including one F2 hybrid between N. 
gabriellae and N. siki. Our results also supported that 
N. leucogenys and N. siki shared the same inversion on 
chromosome 7, which resolves arguments from previous 
studies. Our results demonstrated that both karyotyping 
and DNA-based approaches were suitable for identifying 
purebred species, though neither was ideal for hybrid 
identification. The advantages and disadvantages of 
both approaches are discussed. Our results further 
highlight the importance of animal ethics and welfare, 
which are critical for endangered species in captivity. 
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Keywords: F2 hybrid gibbon; Fluorescence in 
situ hybridization; Nomascus; Pericentric inversion; 
Species identification; Animal welfare 


INTRODUCTION 


The rates of chromosomal evolution are an order of magnitude 
higher in Mammalia than in most other classes of vertebrates 
and are strongly correlated with high speciation rates (Bush 
et al., 1977). High karyotypic diversity and rapid speciation 
are often observed in taxa characterized by small effective 
population sizes, which includes the gibbons and siamangs of 
the small ape family Hylobatidae (Primate). These animals are 
endemic to southern China, South Asia, and Southeast Asia, 
and have been assigned into four genera (Hoolock, Hylobates, 
Nomascus, and Symphalangus; Brandon-Jones et al., 2004; 
Roos & Geissmann, 2001). While siamangs (Symphalangus) 
possess a unique and large gular sac (throat pouch), gibbons 
share very similar external morphology, and were long 
acknowledged as a single genus (Hylobates). | However, 
Hoolock and Nomascus were subsequently recognized as full 
genera based on karyotypic studies, which revealed a unique 
diploid number (2n) of chromosomes among each genus: 
Hoolock (2n=38), Hylobates (2n=44), Nomascus (2n-52), 
and Symphalangus (2n-50) (Chiarelli, 1972; Prouty et al., 
1983; Van Tuinen & Ledbetter, 1983). This subdivision, as 
supported by molecular phylogenetic studies, occurred in the 
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Late Miocene approximately 7 Ma (Fan et al., 2017; Springer 
et al., 2012). Today, 20 species of gibbons and siamangs are 
recognized, though the number continues to increase due to 
newly discovered species (Fan et al., 2017). 

Speciation in this family is considerably faster than that 
observed in many other groups of mammals and is considered 
to be the result of high frequency of chromosomal rearrangement 
(Misceo et al., 2008), including compound inversion/translocation 
(Jauch et al., 1992; Koehler et al., 1995a, 1995b; Marks, 
1982; Nie et al, 2001; Stanyon & Chiarelli, 1982, 1983; 
Van Tuinen & Ledbetter, 1983; Yu et al, 1997). Previous 
comparison between the gibbon (N. leucogenys) and human 
genomes revealed 100 synteny breakpoints, which can facilitate 
chromosome breakage and rearrangement (Carbone et al., 2006). 
Gibbons exhibit highly structured society in their populations, 
which may promote chromosome rearrangement fixation, thereby 
allowing rapid speciation (Wilson et al., 1975). 

The crested gibbon genus Nomascus comprises seven 
recognized species, including ^N. concolor (western 
black-crested gibbon), N.  gabriellae  (yellow-cheeked 
gibbon), N. hainanus (Hainan gibbon), N. leucogenys 
(northern white-cheeked gibbon), N. nasutus (eastern 
black-crested gibbon), N. siki (southern white-cheeked 
gibbon) (Brandon-Jones et al., 2004; Roos & Geissmann, 
2001), and the recently discovered N. annamensis (northern 
buffed-cheeked gibbon) (Thinh et al., 2010). These species 
are characterized by a variety of morphological, anatomical, 
karyological, and vocal features (Couturier & Lernould, 1991; 
Garza & Woodruff, 1994; Schilling, 1984; Thinh et al., 2011). 

Although all studied Nomascus species have a unique 
diploid number of 52 (Couturier et al., 1982; Dutrillaux 
et al., 1975; Van Tuinen & Ledbetter, 1983), they also 
exhibit considerable interspecific karyotypic variation. Using 
the R-banding technique, Couturier & Lernould (1991) 
distinguished the karyotypes for six individuals of WN. 
leucogenys, two N. siki, three N. gabriellae, and several hybrids 
using a combination of an inversion on chromosome 7 and a 
translocation between chromosomes 1 and 22. The inversion 
on chromosome 7 was also confirmed based on the number of 
hybridization signals provided by human chromosome probes 
4 and 22. On unrearranged chromosome 7, one hybridization 
signal of probe 4 on 7q and one signal of probe 22 on 7p were 
detected. On rearranged chromosome 7, two signals of each 
probe were detected, one on each arm of the chromosomes 
(Koehler et al., 1995b). Heterozygous translocation t (1; 22) 
has also been identified in hybrids of N. leucogenys and N. 
gabriellae (Couturier et al., 1982; Turleau et al., 1983). Using 
a PCR-based approach, Carbone et al. (2009) detected 
the inversion on chromosome 7 in N. leucogenys but not in 
N. siki. These results disagreed with those of Couturier & 
Lernould (1991), and thus need to be revisited and clarified 
using additional samples. 

Currently, most gibbon species are endangered or critically 
endangered. Several species, including Hylobates lar and 
N. leucogenys, may have been extirpated from China largely 
due to poaching. Reintroducing animals from zoos to 


their original habitat is a goal for conservation but requires 
accurate identification of purebred animals. Extreme caution 
must be undertaken in captive breeding as interspecific and 
intergeneric hybridization can occur between gibbons and 
siamangs (Hirai et al., 2007; Myers & Shafer, 1979), which 
do not resulting in reduced fertility among offspring in certain 
cases (Van Tuinen et al., 1999). Hybrid offspring are difficult to 
visually identify due to their similar external morphology, with 
identification using the karyotypic approach considered more 
reliable. Co-housing of different species has been avoided in 
many countries and animals are karyotyped for identification, 
with studbooks maintained as a database for appropriate 
management. However, this system is not yet well established 
in China and different species are usually housed together 
resulting in hybridization and unknown interbreeding history. 

In the current study, we examined 21 gibbons using a 
combination of G-banding karyotyping and fluorescence in situ 
hybridization (FISH). We also sequenced mitochondrial and 
nuclear fragments across a known inversion breakpoint on 
chromosome 7 identified in a previous study (Carbone et al., 
2009). We examined whether the combination of G-banding 
and FISH is a reliable approach for species identification. 
We evaluated the efficiency and accuracy of karyotyping- 
and PCR-based approaches for the identification of purebred 
animals and hybrids. We also re-examined whether the 
pericentric inversion on chromosome 7 is species-specific in 
N. leucogenys. 


MATERIALS AND METHODS 


Samples 

We collected peripheral blood from 20 gibbon individuals raised 
in Nanning Zoo, Ningbo Zoo, Kunming Zoo, and Huanglianshan 
National Natural Reserve, and obtained one sample (No.14 in 
Table 1) of cultured gibbon cells maintained at the Kunming 
Institute of Zoology, Chinese Academy of Sciences. Fourteen 
individuals were female and seven were male. All experimental 
procedures and animal care were performed according to the 
protocols approved by the Ethics Committee of the Kunming 
Institute of Zoology, Chinese Academy of Sciences. 


Cell culture, metaphase preparation, and G-banding 
Chromosome suspensions were obtained from lymphocyte 
cultures. Cell culture and metaphase preparations were 
performed following conventional procedures (Hungerford, 
1965). Briefly, whole blood was cultivated in the presence of 
phytohaemagglutinin in RPMI 1640 medium supplemented with 
10% fetal bovine serum. After 68—70 h of growth, colchicine 
was added to the cell cultures to a final concentration of 
0.4—0.8 ug/mL. The cell cultures were incubated for another 
2—4 h before harvest. We treated the cells with hypotonic 
solution (0.075 mol/L KCl) for 20 min, with thrice fixation in 3:1 
methanol/glacial acetic acid. Slides were prepared by applying 
10 uL of metaphase suspension onto dry and clean slides, 
which were then allowed to air dry. G-banding was performed 
following Seabright (1971). All karyotypes were analyzed after 
G-banding. 
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FISH, image capture, and processing 

We used a biotin-labeled probe of the human 22 chromosome. 
Fluorescence in situ hybridization, detection, image capture, 
and processing were carried out following Yang et al. 
(2000) and Nie et al. (2002). We detected fluorescence 
signals using a layer of Cy3-avidin (1:1 000 dilution; 
Amersham Pharmacia Biotech, USA). After detection, slides 


were mounted in Vectashield mounting medium with DAPI 
(4’6-diamidino-2-phenylindole, Vector Laboratories, USA). 
Digital images were acquired using a CytoVision system 
(Applied Imaging Co., USA) with a CCD camera mounted 
on a Zeiss microscope (Germany). We associated the 
hybridization signals with specific chromosome regions based 
on DAPI-banding patterns. 


Table 1 Samples used in this study and a summary of karyotypic characteristics including the FISH signals on the chromosome 7, 


the lengths of chromosomes 1 and 22, the results of PCR using primer sets BP and HSA, as well as the identification results 


based on each approach 





FISH 
Nick Name House (signal on chr. 7) Chr. 1 Chr. 22 Karyotype BP HSA Cyt b genes Final identification 
(pairs) 
L h 
Su-Su Nanning Zoo 2 ONE N. leucogenys N/A N/A N/A N/A 
(Normal) (Normal) 
res Long Short 
Bei-Li Nanning Zoo 2 N. leucogenys * — N. leucogenys N. leucogenys 
(Normal) (Normal) 
Long Short 
Gou-Dan Nanning Zoo 2 N. leucogenys * — N. leucogenys N. leucogenys 
(Normal) (Normal) 
Long Short 
San-Mei Nanning Zoo 2 N. leucogenys N/A N/A N/A N/A 
(Normal) (Normal) 
e Long Short 
No14 Missing 2 N. leucogenys * — N. cf. leucogenys N. cf. leucogenys 
(Normal) (Normal) 
, Long Short 
No name Kunming Zoo 2 N. leucogenys N/A N/A N/A N/A 
(Normal) (Normal) 
, Long Short 
NB1 Ningbo Zoo 2 N. leucogenys N/A N/A N/A N/A 
(Normal) (Normal) 
; Long Short 
NB2 Ningbo Zoo 2 N. leucogenys N/A N/A N/A N/A 
(Normal) (Normal) 
Huanglianshan 
P Long Short 
HLS3 Nature 2 N. leucogenys N/A N/A N/A N/A 
(Normal) (Normal) 
Reserve 
Huanglianshan 
> Long Short 
HLS4 Nature 2 N. leucogenys N/A N/A N/A N/A 
(Normal) (Normal) 
Reserve 
Fang-Fang Nanning Zoo 2 Short Long N. siki N/A N/A N/A N/A 
E’gui Nanning Zoo 2 Short Long N. siki + — N. siki N. siki 
Qingguangyan Nanning Zoo 2 Short Long N. siki N/A N/A N/A N/A 
Laoer Nanning Zoo 2 Short Long N. siki N/A N/A N/A N/A 
N. siki 3 x 
Da-Shan Nanning Zoo 2 Short Long N. siki + — N. gabriellae (N. siki 3 
xN. gabriellae © ) 2 
316 Nanning Zoo 1 Short Long N. gabriellae — N. gabriellae N. gabriellae 
Bai-Shou Nanning Zoo 1 Short Long N. gabriellae — N. gabriellae N. gabriellae 
Jing-Jing Nanning Zoo 1 Short Long N. gabriellae — N. gabriellae N. gabriellae 
1 sh 1 sh N. I 
Mei-Mei Nanning Zoo 2 SrO oR N. leucogenysx N. siki + — N. leucogenys dosi eal 
1 long 1 Long N. siki 
1 N. gabriell. 
Xiao-Xiao Nanning Zoo ane Short Long N. gabriellae x N. siki + + N. leucogenys ii de gen 
a half N. siki 
1 
A-Huang Nanning Zoo z Short Long N. gabriellaex N. siki N/A N/A N/A N/A 
aha 


The karyotype and mitochondrial of Da-Shan suggest different specific affinities. +: The positive result of PCR using primer set BP or HAS. —: 


The negative result of PCR using primer set BP or HAS. N/A: Not available. 
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Amplification, sequencing, and species identification 


After karyotyping, we extracted total DNA from the cultured cells 
using a commercial DNA extraction kit (Blood & Cell Culture 
DNA Mini Kit, Qiagen, Germany). The DNA extracted from 10 
samples was considered acceptable for subsequent analyses as 
other samples were not well preserved after karyotyping. We 
amplified and sequenced cyt b genes using the primer pair 
L14724 hk3 (5/-GGACTTATGACATGAAAAATCATCGTTG-3’) 


and H15915 hk3 (5’-GATTCCCCATTTCTGGTTTACAAGAC-3’). 


We also conducted PCR using the primers provided in Carbone 
et al. (2009) (263C9 BP L with 263C9 BP R (BP hereafter) and 
263C9 BP L with 263C9 HSA22 R (HSA hereafter)), targeting 
a fragment across a breakpoint on chromosome 7. According 
to Carbone et al. (2009), HSA primers should amplify a 
fragment across a breakpoint in species without an inversion 
(e.g., Hylobates spp., N. gabriellae, and N. siki), and BP primers 
should result in positive amplification only for N. leucogenys. The 
mitochondrial and nuclear amplicons were sequenced using a 
ABI 3730 sequencer (Applied Biosystems, USA). 

We identified species based on the karyotyping and 
sequencing results using the following strategy: (1) we 
identified the species based on its karyotype; (2) we verified 
the karyotyping results using PCR with primer sets BP and 
HSA; (3) we BLAST each cyt b gene against the GenBank 
nucleotide database; (4) we estimated a gene tree to verify 
the BLAST results using our cyt b genes accompanied by 
a set of sequences representing the recognized Nomascus 
species downloaded from GenBank; and (5) we repeated 
the karyotyping and amplification/sequencing in triplicate in 
cases where the mitochondrial, nuclear sequencing, and 
karyotyping results conflicted (e.g., sample Da-Shan, see 
Results and Discussion). A final purebred or hybrid and 
species identification was determined. 

We constructed the maximum likelihood gene tree 
using RAxML v8.2.10 and the CIPRES Science Gateway 
(Stamatakis, 2014). In addition to the obtained cyt b 
sequences, 17 cyt b sequences representing seven recognized 
species were downloaded from GenBank (Supplementary 
Table S1) and aligned with our sequences using MAFFT v7.3 
(Katoh & Standley, 2014). We partitioned the alignments 
by codon positions (three partitions). We used GTR+G as 
the evolutionary model for each partition as RAxML does 
not accept models other than GTR or GTR+G. We ran each 
analysis using the rapid bootstrapping algorithm and let RAxML 
halt bootstrapping automatically. 


RESULTS 


Identified karyotypic species 


The diploid number (2n) of all analyzed gibbons was 52. 


G-banding revealed different lengths of chromosomes 1 and 
22, and FISH revealed two, three, or four fluorescence signals 
on chromosome 7 (Figure 1). Considering the G-banded 
karyotyping and FISH results, we recognized five karyotypes 
representing three karyotypic species and two hybrids, with the 
latter two distinguished by nonhomologous pairing. 
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and FISH results with human 22 chromosome-specific probe 
in different Nomascus species and their hybrids 
A: Nomascus leucogenys; B: Nomascus siki, C: Nomascus gabriellae; D: 


hybrid of N. leucogenysxN. siki; E: hybrid of N. gabriellaexN. siki. 


Nomascus leucogenys (Figure 2) 

The G-banded karyotypes of 10 individuals (72, 33) were similar 
to the G-banded (Van Tuinen & Ledbetter, 1983) and R-banded 
karyotypes for N. leucogenys (Couturier & Lernould, 1991). 
Chromosomes 1 and 22 were normal. Two pairs of FISH signals 
were found on chromosome 7, supporting a pericentric inversion 
event (Figure 1A). 
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Figure 2 G-banded karyotype of a northern white-cheeked 
gibbon (N. leucogenys) 


Nomascus siki (Figure 3) 


The karyotypes of five samples (39, 23) were identical to 
those of N. siki reported by Couturier & Lernould (1991). 
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Compared with N. /eucogenys, a balanced translocation 
between chromosome pairs 1 and 22 was detected, resulting 
in one pair of shortened chromosome 1 and one pair of derived 
chromosome 22 (t (1; 22)). Similar to that observed in N. 
leucogenys, FISH demonstrated a pericentric inversion event 
on chromosome 7 (Figure 1B). The sample “Da-Shan” had a 
typical siki-karyotype. 
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Figure 3 G-banded karyotype of a southern white-cheeked 
gibbon (N. siki) 


Nomascus gabriellae (Figure 4) 


The G-banded karyotypes of three animals (1°, 24) were the 
same as that of N. gabriellae (Couturier & Lernould, 1991). 
This karyotype was similar to N. siki in the presence of t (1; 22). 
Only one pair of FISH signals were observed on chromosome 
7, thus differing from N. siki and N. leucogenys (Figure 1C). 
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Figure 4 G-banded karyotype of a yellow-cheeked gibbon (N. 
gabriellae) 


Nomascus leucogenys x Nomascus siki hybrid (Figure 5) 

The G-banded karyotype of one specimen was unique, 
characterized by a heterozygous translocation. One 
chromosome 1 and 22 were normal, similar to that of N. 
leucogenys, and the other chromosome 1 and 22 were similar 
to that of N. siki, indicating one translocation t (1; 22). Two pairs 
of FISH signals were observed on chromosome 7, indicating 
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a pericentric inversion event (Figure 1D). This specimen was 
identified as a hybrid of N. leucogenysx N. siki. 
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Figure 5 G-banded karyotype of a N. leucogenysxN. siki 
hybrid 
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Nomascus gabriellaex Nomascus siki hybrid (Figure 6) 
The distinct G-banded karyotype of two individuals was similar 
to that of N. siki and N. gabriellae in the presence of 
t (1; 22. FISH detected three fluorescence signals on 
chromosome 7, indicating heterozygous pericentric inversion 
(Figure 1E). These two individuals were identified as hybrids of 
N. gabriellaex N. siki. 
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Figure 6 G-banded karyotype of a N. leucogenysxN. siki 
hybrid 


Mitochondrial and nuclear sequences 

We determined the complete cyt b sequences for 10 individuals 
representing the three purebred species (n-7) and two 
hybrids (n=3) identified in the karyotypic analyses. All new 
sequences are available in GenBank under accession Nos. 
MH188408-MH188428 (Supplementary Table S1). Five out 
of the seven purebred animals were supported by BLAST and 
phylogenetic analyses using the cyt b gene and the other two 
were not (Figure 7). One sample (Da-Shan) was identified 
as N. siki based on its karyotype, but its mitochondrial gene 
was typical of N. gabriellae, which also supported by our 
phylogenetic analyses (bootstrap value (BS)=98). Sample 


14 was identified as N. leucogenys based on its karyotype 
but could not be identified unambiguously as it was equally 
similar to N. leucogenys and N. siki. It was closely related 
to (BS=84) and formed a sister lineage of these two species 
on the phylogenetic tree (BS-52). The cyt b genes of the 
two hybrids identified in the karyotypic analyses (Xiao-Xiao 
and Mei-Mei) corresponded to one of their parental species 
(BS>88). On our cyt b genetic tree (Figure 7), sample 14 was 
a close relative to (BS=84) and formed a sister lineage with N. 
siki+N. leucogenys (BS-52). 
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Figure 7 Maximum-likelihood gene tree estimated using the 
cyt b genes of Nomascus species 
Branch lengths correspond to substitution rates. Numbers on branches are 
bootstrap values. Sequences downloaded from GenBank are in gray and 


GenBank accession Nos. are shown. 


The PCR results using the HSA and BP primer sets were 
concordant with the karyotypic results. The PCR results using the 
BP primers (confirming the inversion) were positive for purebred 
N. leucogenys, N. siki and the N. leucogenysx N. siki hybrid 
(Mei-Mei) and were negative when using the HSA primers for 
these samples (Table 1). The PCR results using the BP primers 
were negative for purebred N. gabriellae but positive when using 
the HSA primers. For the hybrid identified as N. gabriellaex N. 
siki (Xiao-Xiao), PCR and sequencing were successful using 
both BP and HSA primers. For Da-Shan, the PCR results were 
positive using BP primers and negative using HSA primers, which 
were congruent with its karyotype (siki type) but disagreed with 


its mitochondrial gene results. 

Because the karyotype and PCR results were congruent, we 
assigned seven samples as purebred species and identified 
another two as hybrids. One sample (14) was identified 
as N. cf. leucogenys as it did not cluster with the other 
N. leucogenys sequences downloaded from GenBank. We 
considered Da-Shan as an F2 hybrid of N. gabriellae and N. 
siki due to the inconsistent patterns revealed by karyotyping 
and PCR using the HSA and BP primers (siki), and the 
mitochondrial phylogeny (gabriellae) (Table 1). 


DISCUSSION 


We used karyotypic- and PCR-based approaches to examine 
the affinities of Nomascus species characterized by high 
karyotypic diversity. We determined the specific karyotypes 
of each species as well as the existence of an inversion on 
chromosome 7 in N. siki and N. leucogenys. We also revealed 
the definitive existence of hybrids in zoos in China, which calls 
attention to the animal ethics and welfare issues related to 
endangered species breeding in captivity. Using the Nomascus 
species as our focal taxa, we evaluated the accuracy and 
efficiency of both approaches in identifying purebred and hybrid 
species. 

Our results supported that N. gabriellae, N. leucogenys, and 
N. siki are distinguishable based on translocation t (1; 22) 
and inversion on chromosome 7 (Figures 1—5), even though N. 
leucogenys and N. siki diverged from each other only recently. 
Couturier & Lernould (1991) revealed that N. hainanus does not 
have an inversion on chromosome 7 or translocation t (1; 22), 
and therefore has a distinct karyotype. It would be interesting 
to examine the G-banded karyotypes and FISH results of 
N. concolor as well as N. annamensis, which is sister to N. 
gabriellae and recognized only recently (Thinh et al., 2010). 
Due to interspecific karyotypic variation, hybrids can be easily 
distinguished based on heterozygous R-, G-banded karyotypes 
and FISH results (Couturier & Lernould, 1991; this study). In 
wild populations, hybridization of Hylobates species such as H. 
larx H. pileatus, H. larx H. agilis, and H. muellerix H. agilis/H. 
albibarbis has been observed (Brockelman & Gittins, 1984). 
Hybridization between Nomascus species in the wild has not 
been reported, which is likely due to their allopatric distribution. 
However, we observed heterozygous karyotypes in captive 
Nomascus animals, congruent with previous study (Couturier 
& Lernould, 1991). Our results support the possibility of 
nonhomologous pairing during meiophase (Figures 5,6), which 
is an interesting characteristic in Hylobatidae. Because 
gibbons have very strong social structures and usually live 
in small family groups, such nonhomologous pairing may be 
easily fixed during evolutionary histories and further promote 
diversification and speciation, resulting in high species diversity 
within the family. 

Our FISH and DNA analyses supported the existence of 
an inversion on chromosome 7 in both N. leucogenys and N. 
siki (Figures 1—3), consistent with the findings of Couturier 
& Lernould (1991), though not supported by Carbone et al. 
(2009). The negative results obtained in the latter study may 
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be due to the hair sample used to represent N. siki, which is 
known for low DNA quantity as well as the existence of PCR 
inhibitors. 

The phylogenetic position of sample 14 was an interesting 
finding. Its karyotype was identical to the other N. leucogenys 
samples (Table 1), which was not supported by the 
mitochondrial gene tree (Figure 7). We repeated amplification 
and sequencing three times for this sample. The sequences 
were identical, indicating no cross-sample contamination, 
and no premature codon was observed, indicating it was 
not a pseudogene. There are two hypothetical scenarios 
that may explain these results. One is incomplete lineage 
sorting between N. leucogenys and N. siki, resulting in 
non-monophyletic relationships. This situation has never been 
observed because the effective population sizes of gibbon 
species are usually small. Alternatively, this sample might 
represent a distinct and unknown taxon, which is a close 
relative but not identical to N. leucogenys or N. siki. The 
sample was obtained from the Kunming Institute of Zoology, 
Chinese Academy of Sciences, in 1993 and the information on 
the animal was not recorded in detail. Unfortunately, neither 
hypothesis could be supported or rejected in the current study. 

Similar to previous studies, our study supported that 
hybridization has occurred in captivity. However, fertility may 
be less impacted as identification of Da-Shan showed it to be 
a likely F2 hybrid of N. gabriellae and N. siki, despite there 
being no available studbook. This is the first record of a 
cross-back F2 hybrid in Nomascus as most other countries 
have well-established systems in place to prevent interbreeding 
of gibbons in captivity, with many animals also previously 
karyotyped. Continuous interbreeding and production of fertile 
offspring may spread across zoos in China because there are 
only limited numbers of gibbons in zoos, and correct species 
identification of their offspring will become far more difficult 
after several generations and recombination. Our findings call 
for the introduction of a system to prevent gibbon interbreeding 
in captivity and for better welfare and awareness of these 
animals. 

The R-banding technique can distinguish different 
karyotypes of Nomascus species (Couturier & Lernould, 
1991), but is both difficult and time-consuming. Herein, we 
demonstrated that G-banding in combination with FISH is 


a reliable approach for identifying karyotypic characteristics. 


This approach was also appropriate for species identification, 
though it was limited in identifying cross-back hybrid F2 


individuals and did not recognize Da-Shan as a hybrid (Table 1). 


As per Carbone et al. (2009), we agree that PCR and 
sequencing can be applicable for examining chromosome 
inversion and species identification. This approach does not 
require high quality cultured cells or karyotyping techniques 
and is applicable for DNA samples with low quality and/or 
low yield, which certainly include, but are not limited to saliva, 
urine, hair, and feces. The known limitation is that accurate 
identification of hybrids relies on a finely established system 
with known karyotypes, breakpoints, and primers. In our case, 
it easily distinguished N. leucogenys/N. siki from N. gabriellae, 
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but could not easily distinguish N. leucogenys/N. siki from the 
hybrid of N. leucogenysx N. siki. 
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ABSTRACT 


Impacts of Quaternary environmental changes on 
mammal faunas of central Asia remain poorly 
understood. due to a lack of comprehensive 
phylogeographic sampling for most species. To 
help address this knowledge gap, we conducted 
the most extensive molecular analysis to date of 
the long-tailed ground squirrel (Urocitellus undulatus 
Pallas 1778) in Mongolia, a country that comprises 
the southern core of this species' range. Drawing on 
material from recent collaborative field expeditions, 
we genotyped 128 individuals at two mitochondrial 
genes (cytochrome b and cytochrome oxidase I; 
1 797 bp total). Phylogenetic inference supports the 
existence of two deeply divergent infraspecific lineages 
(corresponding to subspecies U. u. undulatus and 
U. u. eversmanni), a result in agreement with previous 
molecular investigations but discordant with patterns 
of range-wide craniometric and external phenotypic 
variation. In the widespread western eversmanni 
lineage, we recovered geographically-associated clades 
from the: (a) Khangai, (b) Mongolian Altai, and 
(c) Govi Altai mountain ranges. Phylogeographic 
structure in U. u. eversmanni is consistent with an 
isolation-by- distance model; however, genetic distances 
are significantly lower than among subspecies, and 
intra-clade relationships are largely unresolved. The 
latter patterns, as well as the relatively higher nucleotide 
polymorphism of populations from the Great Lakes 
Depression of northwestern Mongolia, suggest a history 
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of range shifts into these lowland areas in response 
to Pleistocene glaciation and environmental change, 
followed by upslope movements and mitochondrial 
lineage sorting with Holocene aridification. Our study 
illuminates possible historical mechanisms responsible 
for U. undulatus genetic structure and contributes to 
a framework for ongoing exploration of mammalian 
response to past and present climate change in central 
Asia. 


Keywords: Central Asia; Gobi Desert; Great Lakes 
Depression; Mongolia; Phylogeography 


INTRODUCTION 


Urocitellus undulatus Pallas 1778 is a charismatic, 
medium-bodied ground-dwelling sciurid distributed across 
central Asia, including portions of Siberia, Mongolia, 
northwestern China, and easternmost Kazakhstan and 
Kyrgyzstan (Helgen et al., 2009; KryStufek & Vohralik, 2013; 
Ognev, 1947; Wilson & Reeder, 2005). Although the genus 
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Urocitellus (comprised of 12 species formerly subsumed within 
Spermophilus; Helgen et al., 2009) is distributed across much 
of the Holarctic region, U. undulatus is the only exclusively 
Palearctic species in this clade (Wilson & Reeder, 2005; 
McLean et al., 2016b). To date, various single- and multilocus 
investigations (Tsvirka et al., 2008; Ermakov et al., 2015; 
McLean et al., 2016b; Simonov et al., 2017) have revealed that 
U. undulatus is comprised of two deeply divergent lineages 
recognizable as well-defined subspecies (KryStufek & Vohralik, 
2013) or semi-species (Pavlinov & Lissovsky, 2012). These are 
an eastern lineage (undulatus) in western and central Siberia, 
northern Mongolia, and the Amur region of southeastern Siberia, 
and a western lineage (eversmanni) in western Mongolia, 
northern China, and Kazakhstan. However, our understanding 
of the full diversity and evolutionary and biogeographic history of 
this species across its vast range remains incomplete. 

Although U. undulatus has been the subject of persistent 
morphological and molecular focus over the past three decades 
(Ermakov et al., 2015; Linetskaya & Linetskii, 1989; McLean 
et al., 2016b; Simonov et al., 2017; Tsvirka et al., 2008; 
Vorontsov et al., 1980), more expansive genetic datasets 
are necessary to test existing taxonomic hypotheses and 
illuminate the historical demography and biogeography of 
this species. Unfortunately, however, a lack of spatially 
comprehensive sampling and associated genetic data exists 
for this and many other central Asian vertebrate taxa. This 
data gap precludes identification of the broader abiotic and 
biotic processes acting to shape phylogeographic patterns and 


vertebrate community structure across this expansive region. 


For taxa with relatively high morphological conservatism (such 
as Urocitellus), such datasets are particularly crucial to refine 
our understanding of the true genomic and biogeographic 
histories of lineages. 

The most significant lack of phylogeographic sampling from 
U. undulatus is in Mongolia, a country that nevertheless 
comprises the southern core of this species range. Several 
pressing evolutionary and biogeographical questions hinge on 
improved genetic sampling of Mongolian populations. First, 
although each of the subspecific lineages of U. undulatus 
(undulatus and eversmanni) is documented within the country, 
what are their exact geographic distributions? Second, do 
any populations in Mongolia display patterns of mixed mtDNA 
ancestry and, if so, where are these populations located? 
Third, how have known late Quaternary environmental changes 
shaped phylogeographic structure within the widespread 
western lineage (eversmanni)? Specifically, this lineage 
occupies an environmentally and climatically heterogeneous 
range in Mongolia, including multiple mountain systems 
(Khangai, Mongolian Altai, Govi Altai) that were subject to 
late Pleistocene glaciation, downward expansion of permafrost, 
and other environmental changes. 

In this paper, we present the most comprehensive molecular 
phylogeographic analysis of U. undulatus in Mongolia to 
date. Drawing on material collected during expeditions of 
the Museum of Southwestern Biology (New Mexico, USA) 
from 1999—2016, we genotyped samples from across the entire 


Mongolian range of U. undulatus at two mitochondrial genes 
(cytochrome b (cyt b) and cytochrome oxidase | (COl) 1 
797 bp total). We document population genetic variation and 
structure, and use those data to explore the potential effects 
of known late Pleistocene environmental changes on genetic 
patterns. Our work provides new information on the evolutionary 
and biogeographic history of U. undulatus in western Mongolia 
and lays a foundation for further analyses in this and similarly 
distributed central Asian mammals. 


MATERIALS AND METHODS 


Samples and sequencing 

Specimens used in this study are housed at the University 
of New Mexico Museum of Southwestern Biology (MSB). 
Mongolian samples were collected during joint MSB-National 
University of Mongolia expeditions in 1999 (Tinnin et al., 
2002), 2009-2012 (with University of Kansas and University 
of Nebraska), and 2015-2016 (with Northern Michigan 
University). Cumulative efforts of these expeditions include 
>6 500 cataloged mammal specimens from across major 
Mongolian vegetative and faunal provinces, many of which are 
associated with ecto- and endoparasite specimens archived at 
MSB Division of Parasites or University of Nebraska Manter 
Lab of Parasitology. All field methods followed guidelines 
of institutional animal care and use committees as well 
as the American Society of Mammalogists Guide for Use 
of Wild Mammals in Research (Sikes et al., 2011), and 
were focused on collection of “holistic” mammal specimens 
(e.g., Dunnum & Cook, 2012; McLean et al., 2016a; Yates et al., 
1996). Cumulatively, these materials represent an unparalleled 
resource for establishing Mongolian faunal baselines in an era 
of ongoing global climate and environmental change. 

We selected 128 specimens of U. undulatus for sequencing 
and analysis (Figure 1; Supplementary Table S1; GenBank 
accession Nos. MG883400-MG883654). The dataset included 
119 individuals from 11 different Mongolian aimags (political 
land designations analogous to provinces or states) as well as 
putative representatives of both subspecies (as delineated by 
KryStufek & Vohralik, 2013). The dataset also included nine 
individuals of U. u. undulatus from Sakha Republic in northern 
Siberia. We selected sequences of the Columbian ground 
squirrel (U. columbianus) from GenBank as the outgroup 
for phylogenetic analysis. Frozen tissue samples of all 
U. undulatus individuals (liver, muscle, or dried muscle) were 
subjected to lysis in a solution of 600 uL tissue lysis buffer 
and 12-15 uL reconstituted proteinase K (Omega E.Z.N.A. 
kit; Omega Bio-tek, Inc., USA) for up to 24 hours. Genomic 
DNA was isolated using a standard salt/ethanol extraction 
procedure. To reduce the potential for PCR inhibition, all dried 
muscle samples were processed prior to lysis by removing 
debris, cutting into sub-centimeter sized pieces, and washing 
in 10096 ethanol for 15 min at room temperature, vortexing 
several times; these were then washed in STE buffer under 
refrigeration for 12-16 h. Final extractions were quantified 
flourometrically using a Qubit Broad Range assay kit (Life 
Technologies Corp., USA). 
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Figure 1 Map of Mongolia showing major landscape features and sampling localities of U. undulatus 


Higher elevations are shown in darker colors and major mountain ranges and landscape features are indicated with text. Collection localities of samples used in 


this study indicated by red circles (U. u. eversmanni) and squares (U. u. undulatus). 


We used the primer pair MVZO5/MVZ14 (5'-CGAAGC 
TTGATATGAAAAACCATCGTTG/CTTGATATGAAAAACCATCG 
TTG-3/; Smith & Patton, 1993) to amplify all 1 140 bp of the 
mitochondrial cyt b gene. We used the primer pair HCO2198/ 
LCO1490 (5'-TAAACTTCAGGGTGACCAAAAAATCA/GGTC 
AACAAATCATAAAGATATTGG-3'; Folmer et al, 1994) to 
amplify 657 bp of the mitochondrial CO/ gene. Amplification 
of both mtDNA regions took place in 25 uL reactions, with 


annealing temperatures of 52 °C and 48 °C, respectively. 


Purified PCR products were sequenced using Big Dye 
Terminator 3.1 technology (Applied Biosystems, USA) on an 
ABI 3130 automated DNA sequencer in the Molecular Biology 
Facility in the Department of Biology at University of New 
Mexico. Sequences were manually edited in Sequencher 
v5.3 (Gene Codes Corp., Michigan, USA) and aligned with 
MUSCLE v3.7 (Edgar, 2004) using default settings on the 
CIPRES science gateway (www.phylo.org; Miller et al., 2010). 


We used the R packages pegas (Paradis et al., 2017b) and 
popGenome (Pfeifer et al., 2017) to calculate standard population 
genetic statistics, and to test for signals of population expansion 
based on the Tajima's D statistic. Partial deletion of positions 
with missing data was performed when calculating pairwise 
nucleotide-based metrics (nucleotide diversity and pairwise 
number of nucleotide differences). We inferred phylogeny 
for all samples in a Bayesian framework. First, we used 
PartitionFinder v2.1 (Lanfear et al., 2017) to simultaneously 
infer the best-fit partitioning scheme and models of sequence 
evolution for the concatenated mtDNA matrix, as evaluated 
using the AlCc metric. We conducted Bayesian phylogenetic 
inference in MrBayes v3.2.3 (Ronquist & Huelsenbeck, 2003) 
on CIPRES, using the optimal partitioning scheme inferred 
above. Two independent MCMC analyses composed of 4 
Metropolis-coupled chains each (the default) were used to 
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estimate posterior distributions of tree topologies, running both 
analyses for 10 000 000 generations, sampling every 1 000 
generations, and discarding the first 25% of samples as burn-in. 
Convergence of all parameters was assessed in Tracer v1.6.0 
(Rambaut et al., 2014) by visualizing trace plots and ensuring 
effective sample sizes >200. 

To characterize population structure in U. undulatus, we 
performed an analysis of molecular variance (AMOVA) of all 
samples in the R package poppr (Kamvar et al., 2014). For 
each gene, we included subspecies (undulatus, eversmanni) 
and aimag of origin (12 total in our dataset) as factors. We 
also tested significance of the observed variance patterns with 
4 999 randomizations. Because political boundaries may only 
weakly capture landscape-scale genetic structure, we next 
tested consistency of our data with an isolation-by-distance 
(IBD) model, focusing only on U. u. eversmanni. For that 
test, we used functions in the R packages ape (Paradis et al., 
20172) and raster (Hijmans et al., 2017) to compute pairwise 
genetic (P-values) and geographic (in meters) distances, 
respectively. We calculated correlations between these matrices 
using the mantel function in vegan (Oksanen et al., 2017) and 
assessed statistical significance using 4 999 permutations of 
the geographic distance matrix. Finally, we visualized spatial 
patterns of mtDNA diversity by computing a minimum-spanning 
haplotype network (Bandelt et al., 1999) in PopART (http://popart. 
otago.ac.nz), using just the significantly more variable cyt b gene. 


RESULTS 


The best-fit partitioning scheme for the concatenated 
mtDNA matrix included a different partition for each codon 
position within the cyt b and COI genes (although codon 
position 2 for both genes shared the same partition; 
Table 1). The Tamura-Nei (TrN) substitution model or 


one of its extensions (TrN with equal base frequencies, 
gamma-distributed heterogeneity, and/or invariant sites) was 
preferred for all partitions (Table 1). Phylogenetic inference 
in MrBayes recovered two major clades (U. u. undulatus and 
U. u. eversmanni sensu Kryštufek & Vohralik, 2013) with strong 
support (PP=1; Figure 2). The average uncorrected genetic 
distance (mean+SD) between these clades is 5.84%+0.19 for 
cyt b, but a more modest 2.67%+0.20 for COI. Notably, the 


2_Urocitellus_columbianus 
1- Urocitellus columbianus 


U. u. undulatus 


0.0060 
U. u. eversmanni 


0.52 


0.64 


0.84 


maximum inter-clade distance for COI (3.16%) is concordant 
with Ermakov et al.’s (2015) estimate of 3.5% using this 
same marker but different individuals. For reference, average 
uncorrected distances between all samples of U. undulatus 
and the two U. columbianus outgroups are 8.21%+0.14 and 
4.99%+0.16 for cyt b and COI, respectively. However, we note 
that U. undulatus and U. columbianus may not share a most 
recent common ancestor (McLean et al., 2016b). 
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Figure 2 Majority-rule consensus phylogram of Urocitellus undulatus based on Bayesian inference in MrBayes 


Subspecies and major clades (within U. u. eversmanni only) are indicated by text. The tree is rooted on the split from the Nearctic U. columbianus. Branches with 


posterior probabilities less than 0.9 (i.e., between 0.5 and 0.9) are indicated by arrows. Identical cyt b haplotypes are represented by a maximum of 3 individuals 


in the tree; additional duplicate haplotypes were trimmed for clarity (n232). Asterisks denote individuals from Uvs Aimag. 
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Table 1 Best-fit models of evolution according to the AlCc 
metric for partitions of the concatenated mtDNA matrix 








Model No. sites 
cyt b position 1 TrNEF + G 380 
cyt b position 2 + COI position 2 Trn + | 599 
cyt b position 3 TN + G 380 
COI position 1 TrN 219 
CO! position 3 TrN 219 





Both the genotyping results (Table 2) and phylogenetic 
inference confirmed hypotheses that the U. u. eversmanni 
lineage is widespread across western Mongolia, occurring in far 
northern (Khóvsgól), southern (Govi Altai), and westernmost 
(Bayan-Olgii) aimags. Conversely, the nominal eastern lineage 
(U. u. undulatus) occupies a more restricted range in the 
country; individuals taken from as far east as Khóvsgól Aimag, 
Bulgan Aimag and the southeastern Khangai Mountains 
maintain evolutionary affinity with the western U. u. eversmanni 
lineage (Figure 2). We found no evidence for mtDNA admixture 
(i.e., presence of haplotypes from multiple subspecies) in any 
population surveyed, including those proximate to the apparent 
phylogeographic break between undulatus and eversmanni, 
suggesting a persistent lack of gene flow between these two 
subspecific lineages. 

Within the widespread U. u. eversmanni lineage, three 
geographically-associated subclades are strongly supported 
in the MrBayes topology, from (1) the Khangai Mountains 
and surrounding regions (“Khangai” clade); (2) the Mongolian 
and Govi Altai and surrounding highland regions (“Western” 
clade); and (3) additional ranges of the Govi Altai ("Govi" clade). 
Additional, more narrowly distributed genetic clusters were also 
recovered and include populations from Khóvsgól and Uvs 


Aimags (Figure 2), although we note that some clades exhibit 
incomplete lineage sorting. For example, individuals from Uvs 
Aimag are associated with four distinct genetic clusters or 
subclades in the Bayesian phylogeny (asterisks in Figure 2). 
Finer-scale associations with landscape features such as rivers 
were not evident in our dataset. 

Despite the incomplete geographic sorting of mtDNA 
haplotypes and general lack of fine-scale population patterning, 
there is detectable phylogeographic structure within Mongolian 
U. u. eversmanni. AMOVAs support a role for broad provincial 
classifications in mtDNA variation, with 21.7% and 33.9% of 
molecular variance in cyt b and COI, respectively, attributable 
to aimag of origin after accounting for subspecific variation 
(Table 3). Both values are greater than expected by chance 
(Table 3). A statistically significant correlation was also found 
between matrices describing raw genetic distances and raw 
geographic distances within U. u. eversmanni in cyt b (r=0.58, 
P«0.01) and COI (r?=0.38, P«0.01), thereby supporting an 
isolation-by-distance hypothesis in this subspecies. 


Table 2 Population genetic summary statistics for Urocitellus 
undulatus eversmanni, partitioned by gene 
n H Hd .§ k T 
cyt b (1 140bp) 
110 47 0.91 67 9.60 8.71 x 10° 
COI (657bp) 
111 22 0.88 23 2.16 3.29 x 10° 




















n: Number of samples, H: Number of haplotypes, Hd: Haplotype 
diversity, S: Number of polymorphic sites, k: Average number of 
nucleotide differences, x: Nucleotide diversity. 


Table 3 Results of analysis of molecular variance (AMOVA) for both genes and all samples of Urocitellus undulatus 





Degrees 


Sum of 


Variance relative 




















m Variance 
of freedom squared deviations to expected 
cyt b 
Between subspecies 1 2.84 0.05 (8.9) Greater 0.01 
Among aimags 10 15.66 0.12 (21.7) Greater 0.01 
Within aimags 115 42.88 0.37 (69.3) Less 0.01 
Total 126 61.38 0.54 (100) 
COI 
Between subspecies 1 4.25 0.08 (14.9) Greater 0.01 
Among aimags 10 21.33 0.18 (33.9) Greater <0.01 
Within aimags 116 31.65 0.27 (51.2) Less <0.01 
Total 127 57.23 0.53 (100) 





Nevertheless, we emphasize that divergences among 


U. undulatus subclades are very low. This can be visualized 
in the minimum-spanning haplotype network (Figure 3), and 
is borne out in uncorrected pairwise genetic distances 
computed within the clade of (mean+SD) 0.63%+0.39 for cyt 
b and 0.33%+0.22 for COI. Those distances are roughly an 
order of magnitude lower than between U. u. eversmanni 
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and U. u. undulatus (Ermakov et al., 2015; this study). 
They are also lower than those found within U. u. eversmanni 
populations from the adjacent Altai region of southern Russia 
(Simonov et al., 2017), although the latter study used 
the noncoding and more variable mtDNA control region. 
In addition, although there was more variation within than 
among aimags, AMOVAs suggest that there is a significantly 


lower amount of molecular variance within aimags than 
expected by chance, highlighting the shallow differentiation 
that exists in broad geographic regions. Finally, we 
recovered negative values of Tajima’s D for both genes in 


U. u. eversmanni (cyt b, D-—2.09; COI, D=—1.47) and across 
all samples of U. undulatus (cyt b, D-—1.52; COI, D=-0.82), 
although the result was only significant at the P<0.01 level for 
U. u. eversmanni cyt b (P>0.10 for all others). 





Figure 3 Minimum-spanning haplotype network of all Urocitellus undulatus eversmanni samples used in this study 


Individuals in haplotype network (top) are colored by aimag of origin, with colors corresponding to the map (below). Notches indicate single nucleotide 


substitutions. Note the map does not include the distribution of U. u. undulatus (B.-O.-Bayan-Olgii, Govi=Govi-Altai, Khov.=Khovsgol, Bul.=Bulgan, 


Ark.=Arkhangai, Ovor.=Ovorkhangai, Bayan.=Bayankhongor). 


DISCUSSION 


Relatively little is known about range-wide patterns of genetic 
structure and endemicity in many central Asian mammal 
species. This information gap precludes tests of existing 
taxonomic hypotheses and limits deeper knowledge of how 
past environmental changes across this vast region have 
influenced mammalian diversification. This is particularly true 
for ecomorphologically conservative taxa, such as U. undulatus, 
where molecules and morphology might be expected to give 
cryptic or conflicting historical signals. Indeed, the few 
phylogeographic investigations of non-volant vertebrate taxa 
(e.g., Phrynocephalus lizards, Wang & Fu, 2004; Rhombomys 
gerbils, Ning et al., 2007; Bufo toads, Zhang et al., 2008; 


Meriones gerbils and Allactaga, Dipus jerboas, Liao et al., 
2016) in central Asia to date hint at significant impacts of 
late Pleistocene environmental change on population genetic 
diversity and geographic structure in this region. Our analysis 
of U. undulatus allows us to establish preliminary hypotheses 
of mammalian spatiotemporal response to late Quaternary 
change within Mongolia that can be further tested in this and 
other sympatric species. 


On a range-wide scale, our results support the existing 
systematic hypothesis (KryStufek & Vohralik, 2013; Pavlinov 
& Lissovsky, 2012) that U. undulatus is comprised of two 
deeply divergent lineages (undulatus and eversmanni). The 
strong phylogeographic break between these lineages, which 
stretches from southern Lake Baikal (Russia) through eastern 
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Selenge and Tóv Aimags (Mongolia), contrasts, however, 


with previously published patterns of cranial shape variation. 


Specifically, morphological studies (KryStufek & Vohralík, 2013; 
Linetskaya & Linetskii, 1989; Vorontsov et al., 1980) found 
populations from Yakutia and the Amur region of Russia to 
be highly divergent in cranial shape, while populations from 
the remainder of the range extending across southern Siberia 
and northern Mongolia exhibited a broad longitudinal cline 
in cranial shape. Because body size varies significantly in 
U. undulatus, and cranial morphology is highly allometric in 
Urocitellus ground squirrels (e.g., Robinson & Hoffmann, 1975), 
it seems likely that cranial shape data (especially those based 
on linear measurements) largely reflect body size differences, 
which may or may not be useful for elucidating evolutionary 
structure in this species. More robust tests of species limits 
and phylogeographic hypotheses in U. undulatus, as well as 
of our assertion of a lack of gene flow between undulatus 
and eversmanni, await data from additional and independent 
regions of the nuclear genome. 

Considering the eversmanni lineage which makes up 
the bulk of our sampling, our results provide new insights 
into effects of late Quaternary climate change on historical 
biogeography of this taxon across Mongolia. The record of late 
Pleistocene and Holocene environmental change in this region 
includes extensive plateau and mountain valley glaciation, 
specifically in the Khangai, Mongolian Altai, and Govi Altai 
ranges; extensive downward expansion of permafrost; and 
intermittent formation and draining of lakes at both higher and 
lower elevations (Bóhner & Lehmkuhl, 2005; Grunert et al., 
2000; Lehmkuhl, 1998; Lehmkuhl & Lang, 2001). Montane 
glaciation and the expansion of permafrost should have driven 
downslope range shifts in U. undulatus, as this species prefers 
mesic steppe habitats and requires deeper permafrost levels 
for construction of hibernacula. These range shifts, in turn, 
should have promoted increased mixing of populations across 
lowlands of western Mongolia between 30-12 kya. 

Consistent with that scenario, we found shallow divergence 
between major mtDNA clusters within U. u. eversmanni. 
However, because many of these mitochondrial lineages are 
restricted to mid- and high-elevation steppe habitats (e.g., 
in the Govi Altai) and unlikely to experience high levels of 
gene flow, low divergences are likely to be signatures of past 
(i.e., latest Quaternary) population mixing across lowlands of 
western Mongolia. The shallow but significant geographic 
structure that does exist among geographically and ecologically 
disparate populations of U. u. eversmanni could, in turn, 
have been generated by partial lineage sorting and haplotypic 
divergence following expansion into more favorable areas and 
cessation of population connectivity. 

Zhang et al. (2008) found a similar pattern of reduced 
haplotype and nucleotide diversity in green toads (Bufo viridis) 
inhabiting eastern Central Asia. Their data support a history of 
refugial isolation in montane regions of northwest China and 
eastern Kazakhstan followed by rapid postglacial expansion 
into surrounding basins. Liao et al. (2016) described a similar 
pattern in the jerboa Allactaga sibirica in China. Conversely, 
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our data, including low population structure and negative 
values for Tajima's D in U. u. eversmanni, suggest recent 
upslope range expansions from lowland refugia. Therefore, 
from an elevational perspective, these studies support opposite 
historical scenarios that likely reflect differences in the need 
of amphibians to track water availability versus that of steppe 
mammals. Gür (2013) and Liao et al. (2016) describe 
a third pattern in Anatolian ground squirrels (Spermophilus 
xanthoprymnus) in Turkey and gerbils and jerboas (Meriones 
meridianus and Dipus sagitta) in northern China, suggesting 
that these species expanded their areal, but not elevational, 
distributions during the late Pleistocene in conjunction with 
expansion of cold steppe habitats and deserts. 


If our hypothesis of lowland Pleistocene range shifts is 
correct, the most extensive corridors for gene flow among 
ancient populations of U. u. eversmanni may have been 
in the "Great Lakes Depression" and "Valley of the Govi 
Lakes". Those lowland regions, located between major 
Mongolian mountain ranges, are mostly contiguous with a 
broad longitudinal band of mesic steppe that transverses 
Mongolia and forms a corridor for more arid-adapted taxa 
such as Mongolian gazelle (Procapra gutturosa) and Tolai hare 
(Lepus tolai; Batsaikhan et al., 2014). However, Pleistocene 
environments in this region were likely more mesic than today 
and may have included a mixture of steppe and forest elements 
(Grunert et al., 2000; Bóhner & Lehmkuhl, 2005). Downward 
expansion of mesic floral and faunal elements into the Great 
Lakes Depression during the late Pleistocene would have 
provided suitable habitat for an increasingly mesic-adapted 
suite of vertebrate species such as U. undulatus. 


As a post hoc investigation of this scenario, and to more 
thoroughly parse AMOVA results, we calculated population 
genetic statistics (haplotype and nucleotide diversity) for all 
aimags with at least 15 sampled individuals (Uvs, Bayan-Olgii, 
and Govi Altai). While Uvs Aimag contains lower haplotype 
diversity (0.71) than either Bayan-Olgii or Govi Altai aimags 
(0.91 and 0.80, respectively), it has higher nucleotide diversity 
(6.69x103 vs. 2.37x10? and 4.92x10%, respectively). This 
increased nucleotide polymorphism could have resulted from 
confinement of multiple U. u. eversmanni lineages within 
a Pleistocene refugium spread across the Great Lakes 
Depression and surrounding basins, followed by rapid range 
expansion into favorable montane habitats that became 
increasingly disjunct with Holocene climate changes, yielding 
the phylogeographic and demographic signals we detected. 
Simonov et al. (2017) demonstrated elevated mtDNA diversity 
in U. u. eversmanni from the southern Altai Mountains in 
Russia, and suggested that those populations may also have 
been isolated in lowland glacial refugia. Our data strongly 
support their hypothesis that one of these refugia was in the 
Great Lakes Depression. However, we cannot completely 
rule out refugia elsewhere in northern Mongolia, such as 
in Khóvsgól Aimag, a region proximate to the Great Lakes 
Depression, but from which we were only able to sample six 
individuals from a relatively small area. 

Currently, Pleistocene paleoenvironments of western 


Mongolia are somewhat poorly constrained, preventing more 
precise and detailed links between small mammal historical 


biogeography and past environmental change. Grunert et al. 


(2000) proposed a lowstand (i.e., lowered lake levels due to 
local or regional environmental change) for both Uvs Nuur and 
Bayan Nuur during the Last Glacial Maximum (LGM), which 
would have led to exposure of even more extensive areas in 
the Great Lakes Depression than are available today. While 
these lake lowstands are somewhat counterintuitive given 
the relatively mesic conditions inferred for other mid-latitude 
regions during the LGM, a similar pattern of glacial lowstands 
has been described from lakes in nearby northwestern China 
(Fang, 1991). Conversely, the southern Altai Mountains in 
Russia experienced formation of large glacial lakes during 
the late Pleistocene (e.g., Rudoy, 2002). Understanding the 
extent to which these idiosyncratic landscape-level responses 
interacted with regional-scale environmental variability to 
impact the distribution and demography of U. undulatus will 
require linking all currently available sequence datasets with 
new genetic data in a range-wide phylogeographic framework. 


CONCLUSION 


We analyzed the phylogeography of the long tailed ground 
squirrel (Urocitellus undulatus) across the southern core 
of its large central Asian range. Phylogenetic and 
population genetic inferences based on mtDNA strongly 
support the presence of two major lineages in Mongolia 
(U. u. undulatus and U. u. eversmanni). Within the 
more widespread U. u. eversmanni, we identified statistically 
significant but extremely shallow phylogeographic structure, 
with modern genetic clusters associated with Mongolian 
mountain systems (Khangai Mountains, Mongolian and Govi 
Altai). Together, our analyses support a late Pleistocene 
history of extensive population admixture in U. u. eversmanni, 
possibly across the Great Lakes Depression and contiguous 
lowlands of north-west Mongolia, followed by geologically 
recent diversification in postglacial isolation. In addition to 
providing new geographic context for U. undulatus systematics 
and phylogeography, our study establishes hypotheses of 
distributional and demographic response to past environmental 
change in mesic-adapted central Asian mammal species which 
may be tested using robust, genomic-scale datasets. 
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